Changeset 3621

Timestamp:

04/27/10 07:33:58 (15 years ago)

Author:

Eric.Larour

Message:

Almost there

Location:

issm/trunk/src/c

Files:

• DataSet/DataSet.cpp (modified) (3 diffs)
• DataSet/DataSet.h (modified) (1 diff)
• DataSet/Inputs.cpp (modified) (1 diff)
• DataSet/Inputs.h (modified) (2 diffs)
• DataSet/Parameters.cpp (added)
• DataSet/Parameters.h (added)
• EnumDefinitions/EnumDefinitions.h (modified) (2 diffs)
• Makefile.am (modified) (6 diffs)
• ModelProcessorx/CreateDataSets.cpp (modified) (1 diff)
• ModelProcessorx/CreateParameters.cpp (modified) (3 diffs)
• ModelProcessorx/IoModel.h (modified) (11 diffs)
• io/WriteParams.cpp (modified) (1 diff)
• objects/Beam.cpp (modified) (20 diffs)
• objects/Beam.h (modified) (3 diffs)
• objects/Pengrid.cpp (modified) (14 diffs)
• objects/Pengrid.h (modified) (1 diff)
• objects/Penta.cpp (modified) (48 diffs)
• objects/Penta.h (modified) (4 diffs)
• objects/Sing.cpp (modified) (20 diffs)
• objects/Sing.h (modified) (3 diffs)
• objects/Tria.cpp (modified) (2 diffs)
• objects/Tria.h (modified) (2 diffs)

issm/trunk/src/c/DataSet/DataSet.cpp

@@ -392 +392 @@
 }
 /*}}}*/
-/*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<Object*>::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<Object*>::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<Object*>::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<Object*>::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<Object*>::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<Object*>::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<Object*>::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<Object*>::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 +623 @@
         Load* load=NULL;
         Node* node=NULL;
-        Numpar* numpar=NULL;
 
         for ( object=objects.begin() ; object < objects.end(); object++ ){
@@ -886 +640 @@
                         node->Configure(nodes,vertices);
                 }
-                if((*object)->Enum()==NumparEnum){
-                        numpar=(Numpar*)(*object);
-                        numpar->Configure(parameters);
-                }
-
         }
 

issm/trunk/src/c/DataSet/DataSet.h

@@ -48 +48 @@
                 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);

issm/trunk/src/c/DataSet/Inputs.cpp

@@ -249 +249 @@
 }
 /*}}}*/
+/*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<Object*>::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){

issm/trunk/src/c/DataSet/Inputs.h

@@ -6 +6 @@
 #define _INPUTS_H_
 
+class Input;
 #include "./DataSet.h"
 
@@ -28 +29 @@
                 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);
                 /*}}}*/
 

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -99 +99 @@
         InputEnum,
         TriaVertexInputEnum,
+        SingVertexInputEnum,
+        BeamVertexInputEnum,
+        PentaVertexInputEnum,
         BoolInputEnum,
         IntInputEnum,
@@ -155 +158 @@
         ElementOnSurfaceEnum,
         SurfaceAreaEnum,
+        SurfaceSlopexEnum,
+        SurfaceSlopeyEnum,
         WeightsEnum,
         FitEnum,
         AdjointxEnum,
         AdjointyEnum,
+        AdjointzEnum,
+        CollapseEnum,
+        TemperatureEnum,
         PressureEnum
         /*}}}*/

issm/trunk/src/c/Makefile.am

@@ -65 +65 @@
                                         ./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 +87 @@
                                         ./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 +110 @@
                                         ./DataSet/Inputs.cpp\
                                         ./DataSet/Inputs.h\
+                                        ./DataSet/Parameters.cpp\
+                                        ./DataSet/Parameters.h\
                                         ./shared/shared.h\
                                         ./shared/Alloc/alloc.h\
@@ -477 +481 @@
                                         ./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 +503 @@
                                         ./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 +526 @@
                                         ./DataSet/Inputs.cpp\
                                         ./DataSet/Inputs.h\
+                                        ./DataSet/Parameters.cpp\
+                                        ./DataSet/Parameters.h\
                                         ./shared/shared.h\
                                         ./shared/Threads/issm_threads.h\

issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp

@@ -15 +15 @@
 
 
-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: */

issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp

@@ -32 +32 @@
         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 +38 @@
         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 +104 @@
         param->SetDouble(iomodel->dt);
         parameters->AddObject(param);
+
 
         /*ndt: */

issm/trunk/src/c/ModelProcessorx/IoModel.h

@@ -9 +9 @@
 #include "../io/io.h"
 #include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
 #include "../toolkits/toolkits.h"
 
@@ -198 +199 @@
 
         /*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 +209 @@
         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 +232 @@
 
         /*control:*/
-        void  CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+        void  CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
 
         /*thermal:*/
@@ -237 +238 @@
         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 +244 @@
         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 +250 @@
         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 +256 @@
         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 +262 @@
         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 +268 @@
         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 +274 @@
         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);
 
         

issm/trunk/src/c/io/WriteParams.cpp

@@ -63 +63 @@
         for(i=0;i<nfields;i++){
 
-                param=(Param*)parameters->GetObjectByOffset(i+1);
+                param=(Param*)parameters->GetObjectByOffset(i+1); //skip the numpar object
                 
                 switch(param->GetType()){

issm/trunk/src/c/objects/Beam.cpp

@@ -23 +23 @@
 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 +55 @@
                 this->inputs=new Inputs();
         }
-        return;
+        /*point parameters: */
+        this->parameters=beam_parameters;
 }
 /*}}}*/
@@ -68 +69 @@
         int   beam_matice_id;
         int   beam_matpar_id;
-        int   beam_numpar_id;
         int   beam_node_ids[2];
         double nodeinputs[2];
@@ -79 +79 @@
         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 +85 @@
         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 +111 @@
         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 +120 @@
 Beam::~Beam(){
         delete inputs;
-        return;
+        this->parameters=NULL;
 }
 /*}}}*/
@@ -124 +126 @@
 /*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 +135 @@
         hmatice.configure(materialsin);
         hmatpar.configure(materialsin);
-        hnumpar.configure(parametersin);
+
+        /*point parameters to real dataset: */
+        this->parameters=parametersin;
 
 }
@@ -151 +144 @@
 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 +156 @@
         hmatice.DeepEcho();
         hmatpar.DeepEcho();
-        hnumpar.DeepEcho();
+        printf("   parameters\n");
+        parameters->DeepEcho();
         printf("   inputs\n");
         inputs->DeepEcho();
@@ -187 +181 @@
         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 +201 @@
         hmatice.Echo();
         hmatpar.Echo();
-        hnumpar.Echo();
+        printf("   parameters\n");
+        parameters->Echo();
         printf("   inputs\n");
         inputs->Echo();
@@ -236 +233 @@
         hmatice.Marshall(&marshalled_dataset);
         hmatpar.Marshall(&marshalled_dataset);
-        hnumpar.Marshall(&marshalled_dataset);
 
         /*Marshall inputs: */
@@ -243 +239 @@
         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 +255 @@
                 +hmatice.MarshallSize()
                 +hmatpar.MarshallSize()
-                +hnumpar.MarshallSize()
                 +inputs->MarshallSize()
                 +sizeof(int); //sizeof(int) for enum type
@@ -292 +289 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
         
         /*Get dof list on which we will plug the pressure values: */
@@ -301 +297 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*Get node data: */
@@ -454 +449 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -460 +454 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /*recover doflist: */

issm/trunk/src/c/objects/Beam.h

@@ -13 +13 @@
 #include "../ModelProcessorx/IoModel.h"
 #include "../DataSet/Inputs.h"
+#include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
 #include "./Hook.h"
 
@@ -30 +32 @@
                 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 +40 @@
                 /*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();

issm/trunk/src/c/objects/Pengrid.cpp

@@ -18 +18 @@
 #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 +90 @@
         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 +290 @@
 /*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 +298 @@
 /*}}}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 +381 @@
 /*}}}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 +391 @@
         else if (analysis_type==ThermalAnalysisEnum){
                 
-                PenaltyConstrainThermal(punstable,vinputs,analysis_type,sub_analysis_type);
+                PenaltyConstrainThermal(punstable,analysis_type,sub_analysis_type);
                 
         }
@@ -279 +407 @@
 /*}}}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 +488 @@
 /*}}}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 +511 @@
 /*}}}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 +549 @@
 /*}}}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 +597 @@
 /*}}}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 +626 @@
 /*}}}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 +651 @@
 /*}}}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 +723 @@
 /*}}}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;

issm/trunk/src/c/objects/Pengrid.h

@@ -16 +16 @@
 
                 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);
+                /*}}}*/
 
 };

issm/trunk/src/c/objects/Penta.cpp

@@ -17 +17 @@
 #include "../include/typedefs.h"
 #include "../include/macros.h"
+#include "../DataSet/DataSet.h"
+#include "../DataSet/Inputs.h"
 
 /*Object constructors and destructor*/
@@ -22 +24 @@
 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 +56 @@
                 this->inputs=new Inputs();
         }
+        /*point parameters: */
+        this->parameters=penta_parameters;
 }
 /*}}}*/
@@ -64 +67 @@
         int penta_matice_id;
         int penta_matpar_id;
-        int penta_numpar_id;
         double nodeinputs[6];
         bool collapse;
@@ -77 +79 @@
         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 +142 @@
 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 +157 @@
         hmatice.configure(materialsin);
         hmatpar.configure(materialsin);
-        hnumpar.configure(parametersin);
+
+        /*point parameters to real dataset: */
+        this->parameters=parametersin;
 
 }
@@ -174 +165 @@
 /*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 +186 @@
         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 +207 @@
         hmatice.DeepEcho();
         hmatpar.DeepEcho();
-        hnumpar.DeepEcho();
+        printf("   parameters\n");
+        parameters->DeepEcho();
         printf("   inputs\n");
         inputs->DeepEcho();
@@ -230 +224 @@
         hmatice.Echo();
         hmatpar.Echo();
-        hnumpar.Echo();
+        printf("   parameters\n");
+        parameters->Echo();
         printf("   inputs\n");
         inputs->Echo();
@@ -267 +262 @@
         hmatice.Marshall(&marshalled_dataset);
         hmatpar.Marshall(&marshalled_dataset);
-        hnumpar.Marshall(&marshalled_dataset);
 
         /*Marshall inputs: */
@@ -275 +269 @@
         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 +284 @@
                 +hmatice.MarshallSize()
                 +hmatpar.MarshallSize()
-                +hnumpar.MarshallSize()
                 +inputs->MarshallSize()
                 +sizeof(int); //sizeof(int) for enum type
@@ -304 +298 @@
         Hook* tria_hmatice=NULL;
         Hook* tria_hmatpar=NULL;
-        Hook* tria_hnumpar=NULL;
+        Parameters* tria_parameters=NULL;
         Inputs* tria_inputs=NULL;
 
@@ -314 +308 @@
         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 +490 @@
         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 +505 @@
         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 +549 @@
 
                         /*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 +720 @@
         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 +742 @@
 /*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 +771 @@
 /*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 +873 @@
         double MOUNTAINKEXPONENT=10;
 
+        /*parameters: */
+        double viscosity_overshoot;
+
         /*Collapsed formulation: */
         Tria*  tria=NULL;
@@ -866 +882 @@
         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 +898 @@
         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&viscosity_overshoot,"viscosity_overshoot");
 
         /*If on water, skip stiffness: */
@@ -961 +978 @@
                                   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 +1111 @@
         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 +1136 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
 
@@ -1124 +1142 @@
         rho_ice=matpar->GetRhoIce();
         gravity=matpar->GetG();
+
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&stokesreconditioning,"stokesreconditioning");
 
         /* Get node coordinates and dof list: */
@@ -1173 +1194 @@
                         }
                         for (i=6;i<8;i++){
-                                D[i][i]=-D_scalar*numpar->stokesreconditioning;
+                                D[i][i]=-D_scalar*stokesreconditioning;
                         }
 
@@ -1280 +1301 @@
                         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 +1591 @@
 
         int     dofs[3]={0,1,2};
-        double  dt;
         double  K[2][2]={0.0};
 
@@ -1605 +1625 @@
         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 +1652 @@
         nodes=(Node**)hnodes.deliverp();
         matpar=(Matpar*)hmatpar.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -1642 +1663 @@
         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 +1760 @@
 
                         /*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 +1847 @@
         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 +1868 @@
 /*}}}*/
 /*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 +1895 @@
 /*}}}*/
 /*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 +2149 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
+
+        /*parameters: */
+        double stokesreconditioning;
 
         /*inputs: */
@@ -2122 +2163 @@
         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
 
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&stokesreconditioning,"stokesreconditioning");
+
         /*If on water, skip load: */
         if(onwater)return;
@@ -2129 +2173 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
 
@@ -2200 +2243 @@
                         }
                         for (i=6;i<8;i++){
-                                D[i][i]=-D_scalar*numpar->stokesreconditioning;
+                                D[i][i]=-D_scalar*stokesreconditioning;
                         }
 
@@ -2513 +2556 @@
         int       num_vert_gauss=3;
 
-        double dt;
         double temperature_list[numgrids];
         double temperature;
@@ -2552 +2594 @@
         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 +2612 @@
         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
 
+        /*retrieve some parameters: */
+        this->parameters->FindParam(&dt,"dt");
+
         /*If on water, skip: */
         if(onwater)return;
@@ -2575 +2622 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         /* Get node coordinates and dof list: */
@@ -2588 +2634 @@
         beta=matpar->GetBeta();
         meltingpoint=matpar->GetMeltingPoint();
-        dt=numpar->dt;
 
         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 

issm/trunk/src/c/objects/Penta.h

@@ -21 +21 @@
 #include "../ModelProcessorx/IoModel.h"
 #include "./Node.h"
+#include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
+#include "../DataSet/Inputs.h"
 
 class Penta: public Element{
@@ -30 +33 @@
                 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 +41 @@
                 /*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 +107 @@
                 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);

issm/trunk/src/c/objects/Sing.cpp

@@ -16 +16 @@
 #include "../shared/shared.h"
 #include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
 #include "../DataSet/Inputs.h"
 #include "../include/typedefs.h"
@@ -24 +25 @@
 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 +45 @@
 /*}}}*/
 /*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 +59 @@
                 this->inputs=new Inputs();
         }
-        return;
+        /*point parameters: */
+        this->parameters=sing_parameters;
 }
 /*}}}*/
@@ -66 +68 @@
         int sing_matice_id;
         int sing_matpar_id;
-        int sing_numpar_id;
         
         int sing_g;
@@ -79 +80 @@
         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 +85 @@
         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 +91 @@
         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 +100 @@
 Sing::~Sing(){
         delete inputs;
-        return;
+        this->parameters=NULL;
 }
 /*}}}*/
@@ -104 +106 @@
 /*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 +115 @@
         hmatice.configure(materialsin);
         hmatpar.configure(materialsin);
-        hnumpar.configure(parametersin);
+
+        /*point parameters to real dataset: */
+        this->parameters=parametersin;
 
 }
@@ -131 +124 @@
 Object* Sing::copy() {
         
-        return new Sing(*this); 
+        return new Sing(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs);
 
 }
@@ -143 +136 @@
         hmatice.DeepEcho();
         hmatpar.DeepEcho();
-        hnumpar.DeepEcho();
+        printf("   parameters\n");
+        parameters->DeepEcho();
         printf("   inputs\n");
         inputs->DeepEcho();
@@ -167 +161 @@
         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 +179 @@
         hmatice.Echo();
         hmatpar.Echo();
-        hnumpar.Echo();
+        printf("   parameters\n");
+        parameters->Echo();
         printf("   inputs\n");
         inputs->Echo();
@@ -215 +209 @@
         hmatice.Marshall(&marshalled_dataset);
         hmatpar.Marshall(&marshalled_dataset);
-        hnumpar.Marshall(&marshalled_dataset);
 
         /*Marshall inputs: */
@@ -222 +215 @@
         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 +232 @@
                 +hmatice.MarshallSize()
                 +hmatpar.MarshallSize()
-                +hnumpar.MarshallSize()
-                +inputs->MarshallSize()
                 +sizeof(int); //sizeof(int) for enum type
 }
@@ -369 +363 @@
         Matpar* matpar=NULL;
         Matice* matice=NULL;
-        Numpar* numpar=NULL;
 
         /*recover objects from hooks: */
@@ -375 +368 @@
         matpar=(Matpar*)hmatpar.delivers();
         matice=(Matice*)hmatice.delivers();
-        numpar=(Numpar*)hnumpar.delivers();
 
         inputs->GetParameterValue(&slope[0],SurfaceSlopexEnum);

issm/trunk/src/c/objects/Sing.h

@@ -13 +13 @@
 #include "./Hook.h"
 #include "../DataSet/Inputs.h"
+#include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
 
 class Sing: public Element{
@@ -24 +26 @@
                 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 +34 @@
                 /*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();

issm/trunk/src/c/objects/Tria.cpp

@@ -27 +27 @@
 Tria::Tria(){
         this->inputs=NULL;
-        return;
+        this->parameters=NULL;
 }
 /*}}}*/
@@ -148 +148 @@
         delete inputs;
         this->parameters=NULL;
-        return;
 }
 /*}}}*/

issm/trunk/src/c/objects/Tria.h

@@ -6 +6 @@
 #define _TRIA_H_
 
-class Object;
-class Element;
-class Hook;
-class DataSet;
-class Inputs;
-class Node;
-struct IoModel;
-
 #include "./Object.h"
 #include "./Element.h"
@@ -20 +12 @@
 #include "../ModelProcessorx/IoModel.h"
 #include "../DataSet/DataSet.h"
+#include "../DataSet/Parameters.h"
 #include "../DataSet/Inputs.h"