Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3620) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3621) @@ -392,249 +392,4 @@ } /*}}}*/ -/*FUNCTION DataSet::FindParam(double* pscalar, char* name){{{1*/ -int DataSet::FindParam(double* pscalar, char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - int found=0; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Recover the value of this parameter: */ - param->GetParameterValue(pscalar); - found=1; - break; - } - } - } - return found; -} -/*}}}*/ -/*FUNCTION DataSet::FindParam(int* pinteger,char* name){{{1*/ -int DataSet::FindParam(int* pinteger,char* name){ - - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - int found=0; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Recover the value of this parameter: */ - param->GetParameterValue(pinteger); - found=1; - break; - } - } - } - return found; -} -/*}}}*/ -/*FUNCTION DataSet::FindParam(char** pstring,char* name){{{1*/ -int DataSet::FindParam(char** pstring,char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - int found=0; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Recover the value of this parameter: */ - param->GetParameterValue(pstring); - found=1; - break; - } - } - } - return found; - -} -/*}}}*/ -/*FUNCTION DataSet::FindParam(char*** pstringarray,int* pM,char* name){{{1*/ -int DataSet::FindParam(char*** pstringarray,int* pM,char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - int found=0; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Recover the value of this parameter: */ - param->GetParameterValue(pstringarray); - if(pM)*pM=param->GetM(); - found=1; - break; - } - } - } - return found; - -} -/*}}}*/ -/*FUNCTION DataSet::FindParam(double** pdoublearray,int* pM, int* pN,char* name){{{1*/ -int DataSet::FindParam(double** pdoublearray,int* pM, int* pN,char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - int found=0; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Recover the value of this parameter: */ - param->GetParameterValue(pdoublearray); - if(pM)param->GetM(); - if(pN)param->GetN(); - found=1; - break; - } - } - } - return found; - -} -/*}}}*/ -/*FUNCTION DataSet::FindParam(Vec* pvec,char* name){{{1*/ -int DataSet::FindParam(Vec* pvec,char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - int found=0; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Recover the value of this parameter: */ - param->GetParameterValue(pvec); - found=1; - break; - } - } - } - return found; - -} -/*}}}*/ -/*FUNCTION DataSet::FindParamMat* pmat,char* name){{{1*/ -int DataSet::FindParam(Mat* pmat,char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - int found=0; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Recover the value of this parameter: */ - param->GetParameterValue(pmat); - found=1; - break; - } - } - } - return found; - -} -/*}}}*/ -/*FUNCTION DataSet::FindParamObject{{{1*/ -Object* DataSet::FindParamObject(char* name){ - - /*Go through a dataset, and find a Param* object - *which parameter name is "name" : */ - - vector::iterator object; - Param* param=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - /*Find param type objects: */ - if((*object)->Enum()==ParamEnum){ - - /*Ok, this object is a parameter, recover it and ask which name it has: */ - param=(Param*)(*object); - - if (strcmp(param->GetParameterName(),name)==0){ - /*Ok, this is the one! Return the object: */ - return (*object); - } - } - } - return NULL; -} -/*}}}*/ /*FUNCTION DataSet::FindResult(Vec* presult,char* name){{{1*/ int DataSet::FindResult(Vec* presult,char* name){ @@ -868,5 +623,4 @@ Load* load=NULL; Node* node=NULL; - Numpar* numpar=NULL; for ( object=objects.begin() ; object < objects.end(); object++ ){ @@ -886,9 +640,4 @@ node->Configure(nodes,vertices); } - if((*object)->Enum()==NumparEnum){ - numpar=(Numpar*)(*object); - numpar->Configure(parameters); - } - } Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 3620) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 3621) @@ -48,11 +48,5 @@ int DeleteObject(int id); int Size(); - int FindParam(double* pscalar, char* name); - int FindParam(int* pinteger,char* name); - int FindParam(char** pstring,char* name); - int FindParam(char*** pstringarray,int* pM,char* name); - int FindParam(double** pdoublearray,int* pM,int* pN,char* name); - int FindParam(Vec* pvec,char* name); - int FindParam(Mat* pmat,char* name); + int FindResult(Vec* presult,char* name); Object* FindParamObject(char* name); Index: /issm/trunk/src/c/DataSet/Inputs.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3620) +++ /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3621) @@ -249,4 +249,39 @@ } /*}}}*/ +/*FUNCTION Inputs::GetStrainRateStokes(double* epsilon,double* xyz_list, double* gauss, int xenum, int yenum,int zenum){{{1*/ +void Inputs::GetStrainRate(double* epsilon,double* xyz_list, double* gauss, int xenum, int yenum,int zenum){ + + vector::iterator object; + Input* xinput=NULL; + Input* yinput=NULL; + Input* yinput=NULL; + + /*Go through inputs and find data for xenum: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + xinput=(Input*)(*object); + if (xinput->EnumType()==xenum)break; + } + /*Go through inputs and find data for yenum: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + yinput=(Input*)(*object); + if (yinput->EnumType()==yenum)break; + } + /*Go through inputs and find data for zenum: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + zinput=(Input*)(*object); + if (zinput->EnumType()==zenum)break; + } + + if (!xinput | !yinput | !zinput){ + /*we could not find one input with the correct enum type. No defaults values were provided, + * error out: */ + ISSMERROR("%s%i%s%i%s%i\n"," could not find input with enum type ",xenum," or enum type ",yenum, " or enum type ",zenum); + } + + /*Ok, we have the inputs, call bilinear operator: */ + xinput->GetStrainRateStokes(epsilon,yinput,zinput,xyz_list,gauss); + +} +/*}}}*/ /*FUNCTION Inputs::AddInput{{{1*/ int Inputs::AddInput(Input* in_input){ Index: /issm/trunk/src/c/DataSet/Inputs.h =================================================================== --- /issm/trunk/src/c/DataSet/Inputs.h (revision 3620) +++ /issm/trunk/src/c/DataSet/Inputs.h (revision 3621) @@ -6,4 +6,5 @@ #define _INPUTS_H_ +class Input; #include "./DataSet.h" @@ -28,4 +29,5 @@ void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type); void GetStrainRate(double* epsilon,double* xyz_list, double* gauss, int xenum, int yenum); + void GetStrainRateStokes(double* epsilon,double* xyz_list, double* gauss, int xenum, int yenum,int yenum); /*}}}*/ Index: /issm/trunk/src/c/DataSet/Parameters.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Parameters.cpp (revision 3621) +++ /issm/trunk/src/c/DataSet/Parameters.cpp (revision 3621) @@ -0,0 +1,280 @@ +/*!\file Parameters.c + * \brief: implementation of the Parameters object + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./DataSet.h" +#include "./Parameters.h" + +/*Object constructors and destructor*/ +/*FUNCTION Parameters::Parameters(){{{1*/ +Parameters::Parameters(){ + return; +} +/*}}}*/ +/*FUNCTION Parameters::Parameters(int in_enum){{{1*/ +Parameters::Parameters(int in_enum): DataSet(in_enum){ +{ + //do nothing; + return; +} +/*}}}*/ +/*FUNCTION Parameters::~Parameters(){{{1*/ +Parameters::~Parameters(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION Parameters::FindParam(double* pscalar, char* name){{{1*/ +int Parameters::FindParam(double* pscalar, char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pscalar); + found=1; + break; + } + } + } + return found; +} +/*}}}*/ +/*FUNCTION Parameters::FindParam(int* pinteger,char* name){{{1*/ +int Parameters::FindParam(int* pinteger,char* name){ + + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pinteger); + found=1; + break; + } + } + } + return found; +} +/*}}}*/ +/*FUNCTION Parameters::FindParam(char** pstring,char* name){{{1*/ +int Parameters::FindParam(char** pstring,char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pstring); + found=1; + break; + } + } + } + return found; + +} +/*}}}*/ +/*FUNCTION Parameters::FindParam(char*** pstringarray,int* pM,char* name){{{1*/ +int Parameters::FindParam(char*** pstringarray,int* pM,char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pstringarray); + if(pM)*pM=param->GetM(); + found=1; + break; + } + } + } + return found; + +} +/*}}}*/ +/*FUNCTION Parameters::FindParam(double** pdoublearray,int* pM, int* pN,char* name){{{1*/ +int Parameters::FindParam(double** pdoublearray,int* pM, int* pN,char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pdoublearray); + if(pM)param->GetM(); + if(pN)param->GetN(); + found=1; + break; + } + } + } + return found; + +} +/*}}}*/ +/*FUNCTION Parameters::FindParam(Vec* pvec,char* name){{{1*/ +int Parameters::FindParam(Vec* pvec,char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pvec); + found=1; + break; + } + } + } + return found; + +} +/*}}}*/ +/*FUNCTION Parameters::FindParamMat* pmat,char* name){{{1*/ +int Parameters::FindParam(Mat* pmat,char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pmat); + found=1; + break; + } + } + } + return found; + +} +/*}}}*/ +/*FUNCTION Parameters::FindParamObject{{{1*/ +Object* Parameters::FindParamObject(char* name){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Find param type objects: */ + if((*object)->Enum()==ParamEnum){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if (strcmp(param->GetParameterName(),name)==0){ + /*Ok, this is the one! Return the object: */ + return (*object); + } + } + } + return NULL; +} +/*}}}*/ + +/*Object functions*/ Index: /issm/trunk/src/c/DataSet/Parameters.h =================================================================== --- /issm/trunk/src/c/DataSet/Parameters.h (revision 3621) +++ /issm/trunk/src/c/DataSet/Parameters.h (revision 3621) @@ -0,0 +1,30 @@ +/*! \file Parameters.h + * \brief: header file for inputs object + */ + +#ifndef _PARAMETERS_H_ +#define _PARAMETERS_H_ + +#include "./DataSet.h" + +class Parameters: public DataSet{ + + public: + + /*constructors, destructors: {{{1*/ + Parameters(); + Parameters(int enum_type); + ~Parameters(); + /*}}}*/ + /*numerics: {{{1*/ + int FindParam(double* pscalar, char* name); + int FindParam(int* pinteger,char* name); + int FindParam(char** pstring,char* name); + int FindParam(char*** pstringarray,int* pM,char* name); + int FindParam(double** pdoublearray,int* pM,int* pN,char* name); + int FindParam(Vec* pvec,char* name); + int FindParam(Mat* pmat,char* name); + /*}}}*/ + +}; +#endif /* _PARAMETERS_H */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3620) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3621) @@ -99,4 +99,7 @@ InputEnum, TriaVertexInputEnum, + SingVertexInputEnum, + BeamVertexInputEnum, + PentaVertexInputEnum, BoolInputEnum, IntInputEnum, @@ -155,8 +158,13 @@ ElementOnSurfaceEnum, SurfaceAreaEnum, + SurfaceSlopexEnum, + SurfaceSlopeyEnum, WeightsEnum, FitEnum, AdjointxEnum, AdjointyEnum, + AdjointzEnum, + CollapseEnum, + TemperatureEnum, PressureEnum /*}}}*/ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3620) +++ /issm/trunk/src/c/Makefile.am (revision 3621) @@ -65,4 +65,10 @@ ./objects/TriaVertexInput.h\ ./objects/TriaVertexInput.cpp\ + ./objects/SingVertexInput.h\ + ./objects/SingVertexInput.cpp\ + ./objects/BeamVertexInput.h\ + ./objects/BeamVertexInput.cpp\ + ./objects/PentaVertexInput.h\ + ./objects/PentaVertexInput.cpp\ ./objects/BoolInput.h\ ./objects/BoolInput.cpp\ @@ -81,9 +87,5 @@ ./objects/Matpar.h\ ./objects/Matpar.cpp\ - ./objects/Numpar.h\ - ./objects/Numpar.cpp\ ./objects/Input.h\ - ./objects/Input.cpp\ - ./objects/Einput.h\ ./objects/Spc.cpp\ ./objects/Spc.h\ @@ -108,4 +110,6 @@ ./DataSet/Inputs.cpp\ ./DataSet/Inputs.h\ + ./DataSet/Parameters.cpp\ + ./DataSet/Parameters.h\ ./shared/shared.h\ ./shared/Alloc/alloc.h\ @@ -477,4 +481,10 @@ ./objects/TriaVertexInput.h\ ./objects/TriaVertexInput.cpp\ + ./objects/SingVertexInput.h\ + ./objects/SingVertexInput.cpp\ + ./objects/BeamVertexInput.h\ + ./objects/BeamVertexInput.cpp\ + ./objects/PentaVertexInput.h\ + ./objects/PentaVertexInput.cpp\ ./objects/BoolInput.h\ ./objects/BoolInput.cpp\ @@ -493,7 +503,5 @@ ./objects/Matpar.h\ ./objects/Matpar.cpp\ - ./objects/Numpar.h\ - ./objects/Numpar.cpp\ - ./objects/Einput.h\ + ./objects/Input.h\ ./objects/Spc.cpp\ ./objects/Spc.h\ @@ -518,4 +526,6 @@ ./DataSet/Inputs.cpp\ ./DataSet/Inputs.h\ + ./DataSet/Parameters.cpp\ + ./DataSet/Parameters.h\ ./shared/shared.h\ ./shared/Threads/issm_threads.h\ Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 3620) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 3621) @@ -15,5 +15,5 @@ -void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){ /*create parameters common to all solutions: */ Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3620) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3621) @@ -32,9 +32,4 @@ parameters = new DataSet(ParametersEnum); - //solution parameters: always 1'st, to speed-up its lookup by elements.. - count++; - numpar= new Numpar(count); - parameters->AddObject(numpar); - //outputfilename count++; @@ -43,16 +38,4 @@ parameters->AddObject(param); - //analysis and subanalysis - count++; - param= new Param(count,"analysis_type",DOUBLE); - param->SetDouble(iomodel->analysis_type); - parameters->AddObject(param); - - count++; - param= new Param(count,"sub_analysis_type",DOUBLE); - param->SetDouble(iomodel->sub_analysis_type); - parameters->AddObject(param); - - //dimension 2d or 3d: if (strcmp(iomodel->meshtype,"2d")==0)dim=2; @@ -121,4 +104,5 @@ param->SetDouble(iomodel->dt); parameters->AddObject(param); + /*ndt: */ Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3620) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3621) @@ -9,4 +9,5 @@ #include "../io/io.h" #include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" #include "../toolkits/toolkits.h" @@ -198,6 +199,6 @@ /*Creation of fem datasets: general drivers*/ - void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParameters(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); @@ -208,5 +209,5 @@ void CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersDiagnosticHoriz(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*diagnostic vertical*/ @@ -231,5 +232,5 @@ /*control:*/ - void CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*thermal:*/ @@ -237,5 +238,5 @@ void CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersThermal(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersThermal(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*melting:*/ @@ -243,5 +244,5 @@ void CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersMelting(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersMelting(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*prognostic:*/ @@ -249,5 +250,5 @@ void CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersPrognostic(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*prognostic2:*/ @@ -255,5 +256,5 @@ void CreateConstraintsPrognostic2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersPrognostic2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersPrognostic2(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*balancedthickness:*/ @@ -261,5 +262,5 @@ void CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersBalancedthickness(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*balancedthickness2:*/ @@ -267,5 +268,5 @@ void CreateConstraintsBalancedthickness2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersBalancedthickness2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersBalancedthickness2(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*balancedvelocities:*/ @@ -273,8 +274,8 @@ void CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersBalancedvelocities(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersBalancedvelocities(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*qmu: */ - void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersQmu(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); Index: /issm/trunk/src/c/io/WriteParams.cpp =================================================================== --- /issm/trunk/src/c/io/WriteParams.cpp (revision 3620) +++ /issm/trunk/src/c/io/WriteParams.cpp (revision 3621) @@ -63,5 +63,5 @@ for(i=0;iGetObjectByOffset(i+1); + param=(Param*)parameters->GetObjectByOffset(i+1); //skip the numpar object switch(param->GetType()){ Index: /issm/trunk/src/c/objects/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Beam.cpp (revision 3620) +++ /issm/trunk/src/c/objects/Beam.cpp (revision 3621) @@ -23,26 +23,26 @@ Beam::Beam(){ this->inputs=NULL; + this->parameters=NULL; return; } /*}}}*/ -/*FUNCTION Beam::Beam(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id, ElementProperties* properties){{{1*/ -Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id): +/*FUNCTION Beam::Beam(int id, int* node_ids, int matice_id, int matpar_id){{{1*/ +Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id): hnodes(beam_node_ids,2), hmatice(&beam_matice_id,1), - hmatpar(&beam_matpar_id,1), - hnumpar(&beam_numpar_id,1) + hmatpar(&beam_matpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=beam_id; + this->parameters=NULL; this->inputs=new Inputs(); } /*}}}*/ -/*FUNCTION Beam::Beam(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, ElementProperties* properties){{{1*/ -Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, Inputs* beam_inputs): +/*FUNCTION Beam::Beam(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Parameters* beam_parameters, ElementProperties* properties){{{1*/ +Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Parameters* beam_parameters, Inputs* beam_inputs): hnodes(beam_hnodes), hmatice(beam_hmatice), - hmatpar(beam_hmatpar), - hnumpar(beam_hnumpar) + hmatpar(beam_hmatpar) { @@ -55,5 +55,6 @@ this->inputs=new Inputs(); } - return; + /*point parameters: */ + this->parameters=beam_parameters; } /*}}}*/ @@ -68,5 +69,4 @@ int beam_matice_id; int beam_matpar_id; - int beam_numpar_id; int beam_node_ids[2]; double nodeinputs[2]; @@ -79,5 +79,4 @@ beam_matice_id=index+1; //refers to the corresponding material property card beam_matpar_id=iomodel->numberofvertices2d*(iomodel->numlayers-1)+1;//refers to the corresponding matpar property card - beam_numpar_id=1; beam_node_ids[0]=index+1; beam_node_ids[1]=(int)iomodel->uppernodes[index]; //grid that lays right on top @@ -86,5 +85,4 @@ this->hmatice.Init(&beam_matice_id,1); this->hmatpar.Init(&beam_matpar_id,1); - this->hnumpar.Init(&beam_numpar_id,1); //intialize inputs, and add as many inputs per element as requested: @@ -113,4 +111,8 @@ if (iomodel->gridonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->gridonbed[index])); + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; + + } /*}}}*/ @@ -118,5 +120,5 @@ Beam::~Beam(){ delete inputs; - return; + this->parameters=NULL; } /*}}}*/ @@ -124,19 +126,8 @@ /*Object management*/ /*FUNCTION Beam::Configure{{{1*/ -void Beam::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){ +void Beam::Configure(DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){ int i; - DataSet* loadsin=NULL; - DataSet* nodesin=NULL; - DataSet* materialsin=NULL; - DataSet* parametersin=NULL; - - /*Recover pointers :*/ - loadsin=(DataSet*)ploadsin; - nodesin=(DataSet*)pnodesin; - materialsin=(DataSet*)pmaterialsin; - parametersin=(DataSet*)pparametersin; - /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective * datasets, using internal ids and offsets hidden in hooks: */ @@ -144,5 +135,7 @@ hmatice.configure(materialsin); hmatpar.configure(materialsin); - hnumpar.configure(parametersin); + + /*point parameters to real dataset: */ + this->parameters=parametersin; } @@ -151,5 +144,5 @@ Object* Beam::copy() { - return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); + return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs); } @@ -163,5 +156,6 @@ hmatice.DeepEcho(); hmatpar.DeepEcho(); - hnumpar.DeepEcho(); + printf(" parameters\n"); + parameters->DeepEcho(); printf(" inputs\n"); inputs->DeepEcho(); @@ -187,8 +181,10 @@ hmatice.Demarshall(&marshalled_dataset); hmatpar.Demarshall(&marshalled_dataset); - hnumpar.Demarshall(&marshalled_dataset); /*demarshall inputs: */ inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); + + /*parameters: may not exist even yet, so let Configure handle it: */ + this->parameters=NULL; /*return: */ @@ -205,5 +201,6 @@ hmatice.Echo(); hmatpar.Echo(); - hnumpar.Echo(); + printf(" parameters\n"); + parameters->Echo(); printf(" inputs\n"); inputs->Echo(); @@ -236,5 +233,4 @@ hmatice.Marshall(&marshalled_dataset); hmatpar.Marshall(&marshalled_dataset); - hnumpar.Marshall(&marshalled_dataset); /*Marshall inputs: */ @@ -243,4 +239,6 @@ memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); marshalled_dataset+=marshalled_inputs_size; + + /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/ xfree((void**)&marshalled_inputs); @@ -257,5 +255,4 @@ +hmatice.MarshallSize() +hmatpar.MarshallSize() - +hnumpar.MarshallSize() +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type @@ -292,5 +289,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*Get dof list on which we will plug the pressure values: */ @@ -301,5 +297,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*Get node data: */ @@ -454,5 +449,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -460,5 +454,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*recover doflist: */ Index: /issm/trunk/src/c/objects/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Beam.h (revision 3620) +++ /issm/trunk/src/c/objects/Beam.h (revision 3621) @@ -13,4 +13,6 @@ #include "../ModelProcessorx/IoModel.h" #include "../DataSet/Inputs.h" +#include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" #include "./Hook.h" @@ -30,5 +32,6 @@ Hook hmatice; //hook to 1 matice Hook hmatpar; //hook to 1 matpar - Hook hnumpar; //hook to 1 numpar + + Parameters* parameters; //pointer to solution parameters Inputs* inputs; @@ -37,6 +40,6 @@ /*constructors, destructors: {{{1*/ Beam(); - Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id); - Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, Inputs* beam_inputs); + Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id); + Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Parameters* beam_parameters, Inputs* beam_inputs); Beam(int i, IoModel* iomodel); ~Beam(); Index: /issm/trunk/src/c/objects/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Pengrid.cpp (revision 3620) +++ /issm/trunk/src/c/objects/Pengrid.cpp (revision 3621) @@ -18,12 +18,60 @@ #include "../include/typedefs.h" #include "../include/macros.h" - +#include "../DataSet/DataSet.h" +#include "../DataSet/Inputs.h" + /*Object constructors and destructor*/ /*FUNCTION Pengrid::constructor {{{1*/ Pengrid::Pengrid(){ - return; -} -/*}}}1*/ -/*FUNCTION Pengrid::creation {{{1*/ + this->inputs=NULL; + this->parameters=NULL; + + /*not active, not zigzagging: */ + active=0; + zigzag_counter=0; + +} +/*}}}1*/ +/*FUNCTION Pengrid::Pengrid(int id, int node_ids int matpar_id){{{1*/ +Pengrid::Pengrid(int pengrid_id,int pengrid_node_id, int pengrid_matpar_id): + hnode(pengrid_node_ids,1), + hmatice(&pengrid_matice_id,1), + hmatpar(&pengrid_matpar_id,1) +{ + + /*all the initialization has been done by the initializer, just fill in the id: */ + this->id=pengrid_id; + this->parameters=NULL; + this->inputs=new Inputs(); + + /*not active, not zigzagging: */ + active=0; + zigzag_counter=0; + +} +/*}}}*/ +/*FUNCTION Pengrid::Pengrid(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* pengrid_inputs) {{{1*/ +Pengrid::Pengrid(int pengrid_id,Hook* pengrid_hnode, Hook* pengrid_hmatpar, Parameters* pengrid_parameters, Inputs* pengrid_inputs): + hnode(pengrid_hnode), + hmatpar(pengrid_hmatpar) +{ + + /*all the initialization has been done by the initializer, just fill in the id: */ + this->id=pengrid_id; + if(pengrid_inputs){ + this->inputs=(Inputs*)pengrid_inputs->Copy(); + } + else{ + this->inputs=new Inputs(); + } + /*point parameters: */ + this->parameters=pengrid_parameters; + + /*not active, not zigzagging: */ + active=0; + zigzag_counter=0; + +} +/*}}}*/ Pengrid::Pengrid(int pengrid_id, int pengrid_node_id,int pengrid_mparid, int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate,int pengrid_stabilize_constraints){ @@ -42,9 +90,89 @@ matpar_offset=UNDEF; - zigzag_counter=0; return; } /*}}}1*/ +/*FUNCTION Pengrid::Pengrid(int i, IoModel* iomodel){{{1*/ +Pengrid::Pengrid(int index, IoModel* iomodel){ //i is the element index + + int i,j; + int tria_node_ids[3]; + int tria_matice_id; + int tria_matpar_id; + double nodeinputs[3]; + + /*id: */ + this->id=index+1; + + /*hooks: */ + //go recover node ids, needed to initialize the node hook. + if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum){ + /*Discontinuous Galerkin*/ + tria_node_ids[0]=3*index+1; + tria_node_ids[1]=3*index+2; + tria_node_ids[2]=3*index+3; + } + else{ + /*Continuous Galerkin*/ + for(i=0;i<3;i++){ + tria_node_ids[i]=(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab + } + } + tria_matice_id=index+1; //refers to the corresponding ice material object + tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object + + this->hnodes.Init(tria_node_ids,3); + this->hmatice.Init(&tria_matice_id,1); + this->hmatpar.Init(&tria_matpar_id,1); + + //intialize inputs, and add as many inputs per element as requested: + this->inputs=new Inputs(); + + if (iomodel->thickness) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_node_ids[i]-1]; + this->inputs->AddInput(new PengridVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_node_ids[i]-1]; + this->inputs->AddInput(new PengridVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_node_ids[i]-1]; + this->inputs->AddInput(new PengridVertexInput(BedEnum,nodeinputs)); + } + if (iomodel->drag_coefficient) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_node_ids[i]-1]; + this->inputs->AddInput(new PengridVertexInput(DragCoefficientEnum,nodeinputs)); + + if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index])); + if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index])); + this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); + + } + if (iomodel->melting_rate) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_node_ids[i]-1]; + this->inputs->AddInput(new PengridVertexInput(MeltingRateEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_node_ids[i]-1]; + this->inputs->AddInput(new PengridVertexInput(AccumulationRateEnum,nodeinputs)); + } + if (iomodel->geothermalflux) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_node_ids[i]-1]; + this->inputs->AddInput(new PengridVertexInput(GeothermalFluxEnum,nodeinputs)); + } + + if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); + if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); + if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); + if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); + + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; + + +} +/*}}}*/ /*FUNCTION Pengrid::destructor {{{1*/ Pengrid::~Pengrid(){ @@ -162,5 +290,5 @@ /*FUNCTION Pengrid::CreateKMatrix {{{1*/ -void Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Pengrid::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){ /*No loads applied, do nothing: */ @@ -170,5 +298,5 @@ /*}}}1*/ /*FUNCTION Pengrid::CreatePVector {{{1*/ -void Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +void Pengrid::CreatePVector(Vec pg, int analysis_type,int sub_analysis_type){ /*No loads applied, do nothing: */ @@ -253,5 +381,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyConstrain {{{1*/ -void Pengrid::PenaltyConstrain(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyConstrain(int* punstable,int analysis_type,int sub_analysis_type){ if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){ @@ -263,5 +391,5 @@ else if (analysis_type==ThermalAnalysisEnum){ - PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type); + PenaltyConstrainThermal(punstable,analysis_type,sub_analysis_type); } @@ -279,5 +407,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyConstrainThermal {{{1*/ -void Pengrid::PenaltyConstrainThermal(int* punstable,void* vinputs,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyConstrainThermal(int* punstable,int analysis_type,int sub_analysis_type){ // The penalty is stable if it doesn't change during to successive iterations. @@ -360,18 +488,18 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreateMatrix {{{1*/ -void Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyCreateKMatrix(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){ if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){ - PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type); + PenaltyCreateKMatrixDiagnosticStokes( Kgg,kmax,analysis_type,sub_analysis_type); } else if (analysis_type==ThermalAnalysisEnum){ - PenaltyCreateKMatrixThermal( Kgg,inputs,kmax,analysis_type,sub_analysis_type); + PenaltyCreateKMatrixThermal( Kgg,kmax,analysis_type,sub_analysis_type); } else if (analysis_type==MeltingAnalysisEnum){ - PenaltyCreateKMatrixMelting( Kgg,inputs,kmax,analysis_type,sub_analysis_type); + PenaltyCreateKMatrixMelting( Kgg,kmax,analysis_type,sub_analysis_type); } @@ -383,5 +511,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{1*/ -void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){ const int numgrids=1; @@ -421,5 +549,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1*/ -void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){ @@ -469,5 +597,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1*/ -void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){ int found=0; @@ -498,14 +626,14 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreatePVector {{{1*/ -void Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyCreatePVector(Vec pg,double kmax,int analysis_type,int sub_analysis_type){ if (analysis_type==ThermalAnalysisEnum){ - PenaltyCreatePVectorThermal( pg,inputs,kmax,analysis_type,sub_analysis_type); + PenaltyCreatePVectorThermal( pg,kmax,analysis_type,sub_analysis_type); } else if (analysis_type==MeltingAnalysisEnum){ - PenaltyCreatePVectorMelting( pg,inputs,kmax,analysis_type,sub_analysis_type); + PenaltyCreatePVectorMelting( pg,kmax,analysis_type,sub_analysis_type); } @@ -523,5 +651,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{1*/ -void Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyCreatePVectorMelting(Vec pg, double kmax,int analysis_type,int sub_analysis_type){ const int numgrids=1; @@ -595,5 +723,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreatePVectorThermal {{{1*/ -void Pengrid::PenaltyCreatePVectorThermal(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){ +void Pengrid::PenaltyCreatePVectorThermal(Vec pg, double kmax,int analysis_type,int sub_analysis_type){ const int numgrids=1; Index: /issm/trunk/src/c/objects/Pengrid.h =================================================================== --- /issm/trunk/src/c/objects/Pengrid.h (revision 3620) +++ /issm/trunk/src/c/objects/Pengrid.h (revision 3621) @@ -16,53 +16,57 @@ int id; - int dof; - int active; - double penalty_offset; - int thermal_steadystate; - /*nodes: */ - int node_id; - Node* node; - int node_offset; + Hook hnode; //hook to 1 node + Hook hmatpar; //hook to 1 matpar - int mparid; - Matpar* matpar; - int matpar_offset; - - int stabilize_constraints; + Parameters* parameters; //pointer to solution parameters + Inputs* inputs; + + /*internals: */ + int active; int zigzag_counter; public: + /*FUNCTION constructors, destructors {{{1*/ Pengrid(); - Pengrid(int id, int node_id,int mparid,int dof, int active, double penalty_offset,int thermal_steadystate,int stabilize_constraints); + Pengrid(int pengrid_id,int pengrid_node_id int pengrid_matpar_id); + Pengrid(int pengrid_id,Hook* pengrid_hnode, Hook* pengrid_hmatpar, Parameters* pengrid_parameters, Inputs* pengrid_inputs); + Pengrid(int i, IoModel* iomodel); ~Pengrid(); - + /*}}}*/ + /*FUNCTION object management {{{1*/ + void Configure(void* elements,void* nodes,void* materials); + Object* copy(); + void DeepEcho(); + void Demarshall(char** pmarshalled_dataset); void Echo(); - void DeepEcho(); + int Enum(); + int GetId(); + char* GetName(); void Marshall(char** pmarshalled_dataset); int MarshallSize(); - char* GetName(); - void Demarshall(char** pmarshalled_dataset); - int Enum(); - int GetId(); int MyRank(); + /*}}}*/ + /*FUNCTION element numerical routines {{{1*/ void DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type); - void Configure(void* elements,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); - void UpdateFromInputs(void* inputs); - void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); - void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); - void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); + void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, int analysis_type,int sub_analysis_type); + void PenaltyCreateKMatrix(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVector(Vec pg,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofspernode); - Object* copy(); - void PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type); - void PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreateKMatrixThermal(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreateKMatrixMelting(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type); void MatparConfiguration(Matpar* matpar,int matpar_offset); - void PenaltyCreatePVectorThermal(Vec pg, void* inputs, double kmax,int analysis_type,int sub_analysis_type); - void PenaltyCreatePVectorMelting(Vec pg, void* inputs, double kmax,int analysis_type,int sub_analysis_type); - void PenaltyConstrain(int* punstable,void* inputs,int analysis_type,int sub_analysis_type); - void PenaltyConstrainThermal(int* punstable,void* inputs,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVectorThermal(Vec pg, double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVectorMelting(Vec pg, double kmax,int analysis_type,int sub_analysis_type); + void PenaltyConstrain(int* punstable,int analysis_type,int sub_analysis_type); + void PenaltyConstrainThermal(int* punstable,int analysis_type,int sub_analysis_type); + + /*updates:*/ + void UpdateFromDakota(void* inputs); + void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); + /*}}}*/ }; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 3620) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 3621) @@ -17,4 +17,6 @@ #include "../include/typedefs.h" #include "../include/macros.h" +#include "../DataSet/DataSet.h" +#include "../DataSet/Inputs.h" /*Object constructors and destructor*/ @@ -22,27 +24,26 @@ Penta::Penta(){ this->inputs=NULL; - return; + this->parameters=NULL; } /*}}}*/ /*FUNCTION Penta constructor {{{1*/ -Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id): +Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id): hnodes(penta_node_ids,6), hmatice(&penta_matice_id,1), - hmatpar(&penta_matpar_id,1), - hnumpar(&penta_numpar_id,1) + hmatpar(&penta_matpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=penta_id; + this->parameters=NULL; this->inputs=new Inputs(); } /*}}}*/ -/*FUNCTION Penta other constructor {{{1*/ -Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, Inputs* penta_inputs): +/*FUNCTION Penta::Penta(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* penta_inputs) {{{1*/ +Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Parameters* penta_parameters, Inputs* penta_inputs): hnodes(penta_hnodes), hmatice(penta_hmatice), - hmatpar(penta_hmatpar), - hnumpar(penta_hnumpar) + hmatpar(penta_hmatpar) { @@ -55,4 +56,6 @@ this->inputs=new Inputs(); } + /*point parameters: */ + this->parameters=penta_parameters; } /*}}}*/ @@ -64,5 +67,4 @@ int penta_matice_id; int penta_matpar_id; - int penta_numpar_id; double nodeinputs[6]; bool collapse; @@ -77,10 +79,8 @@ penta_matice_id=index+1; //refers to the corresponding ice material object penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object - penta_numpar_id=1; //refers to numerical parameters object this->hnodes.Init(&penta_node_ids[0],6); this->hmatice.Init(&penta_matice_id,1); this->hmatpar.Init(&penta_matpar_id,1); - this->hnumpar.Init(&penta_numpar_id,1); //intialize inputs, and add as many inputs per element as requested: @@ -142,24 +142,13 @@ Penta::~Penta(){ delete inputs; - return; + this->parameters=NULL; } /*}}}*/ /*Object management: */ -/*FUNCTION Configure {{{1*/ -void Penta::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){ +/*FUNCTION Penta::Configure {{{1*/ +void Penta::Configure(DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){ int i; - - DataSet* loadsin=NULL; - DataSet* nodesin=NULL; - DataSet* materialsin=NULL; - DataSet* parametersin=NULL; - - /*Recover pointers :*/ - loadsin=(DataSet*)ploadsin; - nodesin=(DataSet*)pnodesin; - materialsin=(DataSet*)pmaterialsin; - parametersin=(DataSet*)pparametersin; /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective @@ -168,5 +157,7 @@ hmatice.configure(materialsin); hmatpar.configure(materialsin); - hnumpar.configure(parametersin); + + /*point parameters to real dataset: */ + this->parameters=parametersin; } @@ -174,5 +165,5 @@ /*FUNCTION copy {{{1*/ Object* Penta::copy() { - return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); + return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs); } /*}}}*/ @@ -195,8 +186,10 @@ hmatice.Demarshall(&marshalled_dataset); hmatpar.Demarshall(&marshalled_dataset); - hnumpar.Demarshall(&marshalled_dataset); /*demarshall inputs: */ inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); + + /*parameters: may not exist even yet, so let Configure handle it: */ + this->parameters=NULL; /*return: */ @@ -214,5 +207,6 @@ hmatice.DeepEcho(); hmatpar.DeepEcho(); - hnumpar.DeepEcho(); + printf(" parameters\n"); + parameters->DeepEcho(); printf(" inputs\n"); inputs->DeepEcho(); @@ -230,5 +224,6 @@ hmatice.Echo(); hmatpar.Echo(); - hnumpar.Echo(); + printf(" parameters\n"); + parameters->Echo(); printf(" inputs\n"); inputs->Echo(); @@ -267,5 +262,4 @@ hmatice.Marshall(&marshalled_dataset); hmatpar.Marshall(&marshalled_dataset); - hnumpar.Marshall(&marshalled_dataset); /*Marshall inputs: */ @@ -275,6 +269,7 @@ marshalled_dataset+=marshalled_inputs_size; + /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/ + xfree((void**)&marshalled_inputs); - *pmarshalled_dataset=marshalled_dataset; @@ -289,5 +284,4 @@ +hmatice.MarshallSize() +hmatpar.MarshallSize() - +hnumpar.MarshallSize() +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type @@ -304,5 +298,5 @@ Hook* tria_hmatice=NULL; Hook* tria_hmatpar=NULL; - Hook* tria_hnumpar=NULL; + Parameters* tria_parameters=NULL; Inputs* tria_inputs=NULL; @@ -314,13 +308,12 @@ tria_hmatice=this->hmatice.Spawn(&zero,1); tria_hmatpar=this->hmatpar.Spawn(&zero,1); - tria_hnumpar=this->hnumpar.Spawn(&zero,1); + tria_parameters=this->parameters; tria_inputs=(Inputs*)this->inputs->Spawn(indices,3); - tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_inputs); + tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_parameters,tria_inputs); delete tria_hnodes; delete tria_hmatice; delete tria_hmatpar; - delete tria_hnumpar; delete tria_inputs; @@ -497,9 +490,11 @@ bool onbed; + /*parameters: */ + double stokesreconditioning; + /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*Check analysis_types*/ @@ -510,9 +505,11 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*recover some inputs: */ inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*retrieve some parameters: */ + this->parameters->FindParam(&stokesreconditioning,"stokesreconditioning"); + if(!onbed){ //put zero @@ -552,7 +549,7 @@ /*Compute Stress*/ - sigma_xx=viscosity*epsilon[0]-pressure*numpar->stokesreconditioning; // sigma = nu eps - pressure - sigma_yy=viscosity*epsilon[1]-pressure*numpar->stokesreconditioning; - sigma_zz=viscosity*epsilon[2]-pressure*numpar->stokesreconditioning; + sigma_xx=viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure + sigma_yy=viscosity*epsilon[1]-pressure*stokesreconditioning; + sigma_zz=viscosity*epsilon[2]-pressure*stokesreconditioning; sigma_xy=viscosity*epsilon[3]; sigma_xz=viscosity*epsilon[4]; @@ -723,9 +720,9 @@ else if (analysis_type==BalancedthicknessAnalysisEnum){ - CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixBalancedthickness( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixBalancedvelocities( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==ThermalAnalysisEnum){ @@ -745,18 +742,26 @@ /*FUNCTION CreateKMatrixBalancedthickness {{{1*/ -void Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixBalancedthickness(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ Tria* tria=NULL; + /*flags: */ + bool onwater; + bool onbed; + + /*recover some inputs: */ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type); delete tria; return; @@ -766,18 +771,26 @@ /*FUNCTION CreateKMatrixBalancedvelocities {{{1*/ -void Penta::CreateKMatrixBalancedvelocities(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixBalancedvelocities(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ Tria* tria=NULL; + /*flags: */ + bool onbed; + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrix(Kgg,analysis_type,sub_analysis_type); delete tria; return; @@ -860,4 +873,7 @@ double MOUNTAINKEXPONENT=10; + /*parameters: */ + double viscosity_overshoot; + /*Collapsed formulation: */ Tria* tria=NULL; @@ -866,10 +882,8 @@ Node** nodes=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*inputs: */ @@ -884,4 +898,7 @@ inputs->GetParameterValue(&onbed,ElementOnBedEnum); inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&viscosity_overshoot,"viscosity_overshoot"); /*If on water, skip stiffness: */ @@ -961,5 +978,5 @@ onto this scalar matrix, so that we win some computational time: */ - newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity); + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=newviscosity*gauss_weight*Jdet; for (i=0;i<5;i++){ @@ -1094,9 +1111,11 @@ double drag_p,drag_q; + /*parameters: */ + double stokesreconditioning; + /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*inputs: */ @@ -1117,5 +1136,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); @@ -1124,4 +1142,7 @@ rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); + + /*retrieve some parameters: */ + this->parameters->FindParam(&stokesreconditioning,"stokesreconditioning"); /* Get node coordinates and dof list: */ @@ -1173,5 +1194,5 @@ } for (i=6;i<8;i++){ - D[i][i]=-D_scalar*numpar->stokesreconditioning; + D[i][i]=-D_scalar*stokesreconditioning; } @@ -1280,7 +1301,7 @@ DLStokes[9][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0; DLStokes[10][10]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0; - DLStokes[11][11]=numpar->stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0]; - DLStokes[12][12]=numpar->stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1]; - DLStokes[13][13]=numpar->stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2]; + DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0]; + DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1]; + DLStokes[13][13]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2]; /* Do the triple product tL*D*L: */ @@ -1570,5 +1591,4 @@ int dofs[3]={0,1,2}; - double dt; double K[2][2]={0.0}; @@ -1605,8 +1625,10 @@ double mixed_layer_capacity,thermal_exchange_velocity; + /*parameters: */ + double dt,artdiff,epsvel; + /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; Matpar* matpar=NULL; - Numpar* numpar=NULL; /*Collapsed formulation: */ @@ -1630,5 +1652,4 @@ nodes=(Node**)hnodes.deliverp(); matpar=(Matpar*)hmatpar.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -1642,5 +1663,9 @@ heatcapacity=matpar->GetHeatCapacity(); thermalconductivity=matpar->GetThermalConductivity(); - dt=numpar->dt; + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); + this->parameters->FindParam(&artdiff,"art_diff"); + this->parameters->FindParam(&epsvel,"epsvel"); /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore @@ -1735,7 +1760,7 @@ /*Artifficial diffusivity*/ - if(numpar->artdiff){ + if(artdiff){ /*Build K: */ - D_scalar=gauss_weight*Jdet/(pow(u,2)+pow(v,2)+numpar->epsvel); + D_scalar=gauss_weight*Jdet/(pow(u,2)+pow(v,2)+epsvel); if(dt){ D_scalar=D_scalar*dt; @@ -1822,9 +1847,9 @@ else if (analysis_type==BalancedthicknessAnalysisEnum){ - CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorPrognostic( pg,analysis_type,sub_analysis_type); } else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorPrognostic( pg,analysis_type,sub_analysis_type); } else if (analysis_type==ThermalAnalysisEnum){ @@ -1843,19 +1868,26 @@ /*}}}*/ /*FUNCTION CreatePVectorBalancedthickness {{{1*/ - -void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorBalancedthickness( Vec pg, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ Tria* tria=NULL; + /*flags: */ + bool onbed; + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVector(pg, analysis_type,sub_analysis_type); delete tria; return; @@ -1863,19 +1895,26 @@ /*}}}*/ /*FUNCTION CreatePVectorBalancedvelocities {{{1*/ - -void Penta::CreatePVectorBalancedvelocities( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorBalancedvelocities( Vec pg, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ Tria* tria=NULL; + /*flags: */ + bool onbed; + bool onwater; + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVector(pg, analysis_type,sub_analysis_type); delete tria; return; @@ -2110,5 +2149,7 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; + + /*parameters: */ + double stokesreconditioning; /*inputs: */ @@ -2122,4 +2163,7 @@ inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + /*retrieve some parameters: */ + this->parameters->FindParam(&stokesreconditioning,"stokesreconditioning"); + /*If on water, skip load: */ if(onwater)return; @@ -2129,5 +2173,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); @@ -2200,5 +2243,5 @@ } for (i=6;i<8;i++){ - D[i][i]=-D_scalar*numpar->stokesreconditioning; + D[i][i]=-D_scalar*stokesreconditioning; } @@ -2513,5 +2556,4 @@ int num_vert_gauss=3; - double dt; double temperature_list[numgrids]; double temperature; @@ -2552,9 +2594,11 @@ Tria* tria=NULL; + /*parameters: */ + double dt; + /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*inputs: */ @@ -2568,4 +2612,7 @@ inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); + /*If on water, skip: */ if(onwater)return; @@ -2575,5 +2622,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -2588,5 +2634,4 @@ beta=matpar->GetBeta(); meltingpoint=matpar->GetMeltingPoint(); - dt=numpar->dt; /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 3620) +++ /issm/trunk/src/c/objects/Penta.h (revision 3621) @@ -21,4 +21,7 @@ #include "../ModelProcessorx/IoModel.h" #include "./Node.h" +#include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" +#include "../DataSet/Inputs.h" class Penta: public Element{ @@ -30,6 +33,6 @@ Hook hmatice; //hook to 1 matice Hook hmatpar; //hook to 1 matpar - Hook hnumpar; //hook to 1 numpar + Parameters* parameters; //pointer to solution parameters Inputs* inputs; @@ -38,6 +41,6 @@ /*FUNCTION constructors, destructors {{{1*/ Penta(); - Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id); - Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, Inputs* inputs); + Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id); + Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Parameters* penta_parameters, Inputs* inputs); Penta(int i, IoModel* iomodel); ~Penta(); @@ -104,6 +107,9 @@ void CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type); void CreatePVectorPrognostic( Vec pg, int analysis_type,int sub_analysis_type); - + void CreateKMatrixBalancedthickness(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreateKMatrixBalancedvelocities(Mat Kgg,int analysis_type,int sub_analysis_type); void CreateKMatrixDiagnosticStokes( Mat Kgg, int analysis_type,int sub_analysis_type); + void CreatePVectorBalancedthickness( Vec pg, int analysis_type,int sub_analysis_type); + void CreatePVectorBalancedvelocities( Vec pg, int analysis_type,int sub_analysis_type); void CreatePVectorDiagnosticStokes( Vec pg, int analysis_type,int sub_analysis_type); void ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp); Index: /issm/trunk/src/c/objects/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Sing.cpp (revision 3620) +++ /issm/trunk/src/c/objects/Sing.cpp (revision 3621) @@ -16,4 +16,5 @@ #include "../shared/shared.h" #include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" #include "../DataSet/Inputs.h" #include "../include/typedefs.h" @@ -24,17 +25,18 @@ Sing::Sing(){ this->inputs=NULL; + this->parameters=NULL; return; } /*}}}*/ /*FUNCTION Sing constructor {{{1*/ -Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id): +Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id): hnodes(sing_node_ids,1), hmatice(&sing_matice_id,1), - hmatpar(&sing_matpar_id,1), - hnumpar(&sing_numpar_id,1) + hmatpar(&sing_matpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=sing_id; + this->parameters=NULL; this->inputs=new Inputs(); @@ -43,9 +45,8 @@ /*}}}*/ /*FUNCTION Sing other constructor {{{1*/ -Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs): +Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Parameters* sing_parameters,Inputs* sing_inputs): hnodes(sing_hnodes), hmatice(sing_hmatice), - hmatpar(sing_hmatpar), - hnumpar(sing_hnumpar) + hmatpar(sing_hmatpar) { @@ -58,5 +59,6 @@ this->inputs=new Inputs(); } - return; + /*point parameters: */ + this->parameters=sing_parameters; } /*}}}*/ @@ -66,5 +68,4 @@ int sing_matice_id; int sing_matpar_id; - int sing_numpar_id; int sing_g; @@ -79,5 +80,4 @@ sing_matice_id=i+1; //refers to the corresponding material property card sing_matpar_id=iomodel->numberofvertices+1;//refers to the corresponding matpar property card - sing_numpar_id=1; sing_g=i+1; @@ -85,5 +85,4 @@ this->hmatice.Init(&sing_matice_id,1); this->hmatpar.Init(&sing_matpar_id,1); - this->hnumpar.Init(&sing_numpar_id,1); //intialize inputs, and add as many inputs per element as requested: @@ -92,4 +91,7 @@ if (iomodel->thickness) this->inputs->AddInput(new SingVertexInput(ThicknessEnum,iomodel->thickness[i])); if (iomodel->drag_coefficient) this->inputs->AddInput(new SingVertexInput(DragCoefficientEnum,iomodel->drag_coefficient[i])); + + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; } @@ -98,5 +100,5 @@ Sing::~Sing(){ delete inputs; - return; + this->parameters=NULL; } /*}}}*/ @@ -104,19 +106,8 @@ /*Object management*/ /*FUNCTION Sing::Configure {{{1*/ -void Sing::Configure(void* ploadsin, void* pnodesin,void* pmaterialsin,void* pparametersin){ +void Sing::Configure(DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){ int i; - DataSet* loadsin=NULL; - DataSet* nodesin=NULL; - DataSet* materialsin=NULL; - DataSet* parametersin=NULL; - - /*Recover pointers :*/ - loadsin=(DataSet*)ploadsin; - nodesin=(DataSet*)pnodesin; - materialsin=(DataSet*)pmaterialsin; - parametersin=(DataSet*)pparametersin; - /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective * datasets, using internal ids and offsets hidden in hooks: */ @@ -124,5 +115,7 @@ hmatice.configure(materialsin); hmatpar.configure(materialsin); - hnumpar.configure(parametersin); + + /*point parameters to real dataset: */ + this->parameters=parametersin; } @@ -131,5 +124,5 @@ Object* Sing::copy() { - return new Sing(*this); + return new Sing(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs); } @@ -143,5 +136,6 @@ hmatice.DeepEcho(); hmatpar.DeepEcho(); - hnumpar.DeepEcho(); + printf(" parameters\n"); + parameters->DeepEcho(); printf(" inputs\n"); inputs->DeepEcho(); @@ -167,8 +161,7 @@ hmatice.Demarshall(&marshalled_dataset); hmatpar.Demarshall(&marshalled_dataset); - hnumpar.Demarshall(&marshalled_dataset); - - /*demarshall inputs: */ - inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); + + /*parameters: may not exist even yet, so let Configure handle it: */ + this->parameters=NULL; /*return: */ @@ -186,5 +179,6 @@ hmatice.Echo(); hmatpar.Echo(); - hnumpar.Echo(); + printf(" parameters\n"); + parameters->Echo(); printf(" inputs\n"); inputs->Echo(); @@ -215,5 +209,4 @@ hmatice.Marshall(&marshalled_dataset); hmatpar.Marshall(&marshalled_dataset); - hnumpar.Marshall(&marshalled_dataset); /*Marshall inputs: */ @@ -222,4 +215,7 @@ memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); marshalled_dataset+=marshalled_inputs_size; + + /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/ + xfree((void**)&marshalled_inputs); @@ -236,6 +232,4 @@ +hmatice.MarshallSize() +hmatpar.MarshallSize() - +hnumpar.MarshallSize() - +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type } @@ -369,5 +363,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -375,5 +368,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); inputs->GetParameterValue(&slope[0],SurfaceSlopexEnum); Index: /issm/trunk/src/c/objects/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Sing.h (revision 3620) +++ /issm/trunk/src/c/objects/Sing.h (revision 3621) @@ -13,4 +13,6 @@ #include "./Hook.h" #include "../DataSet/Inputs.h" +#include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" class Sing: public Element{ @@ -24,5 +26,6 @@ Hook hmatice; //hook to 1 matice Hook hmatpar; //hook to 1 matpar - Hook hnumpar; //hook to 1 numpar + + Parameters* parameters; //pointer to solution parameters Inputs* inputs; @@ -31,6 +34,6 @@ /*constructors, destructors: {{{1*/ Sing(); - Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id); - Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs); + Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id); + Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Parameters* sing_parameters,Inputs* sing_inputs); Sing(int i, IoModel* iomodel); ~Sing(); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3620) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3621) @@ -27,5 +27,5 @@ Tria::Tria(){ this->inputs=NULL; - return; + this->parameters=NULL; } /*}}}*/ @@ -148,5 +148,4 @@ delete inputs; this->parameters=NULL; - return; } /*}}}*/ Index: /issm/trunk/src/c/objects/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Tria.h (revision 3620) +++ /issm/trunk/src/c/objects/Tria.h (revision 3621) @@ -6,12 +6,4 @@ #define _TRIA_H_ -class Object; -class Element; -class Hook; -class DataSet; -class Inputs; -class Node; -struct IoModel; - #include "./Object.h" #include "./Element.h" @@ -20,4 +12,5 @@ #include "../ModelProcessorx/IoModel.h" #include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" #include "../DataSet/Inputs.h"