Changeset 246


Ignore:
Timestamp:
05/05/09 14:11:46 (16 years ago)
Author:
Eric.Larour
Message:

New parameter inputs framework, that can adapt to unknow number of dofs problems

Location:
issm/trunk/src
Files:
12 added
45 edited

Legend:

Unmodified
Added
Removed

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

    r202 r246  
    881881#undef __FUNCT__
    882882#define __FUNCT__ "DataSet::CreateKMatrix"
    883 void  DataSet::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     883void  DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
    884884
    885885        vector<Object*>::iterator object;
    … …  
    905905#undef __FUNCT__
    906906#define __FUNCT__ "DataSet::CreatePVector"
    907 void  DataSet::CreatePVector(Vec pg,ParameterInputs* inputs,int analysis_type){
     907void  DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type){
    908908
    909909        vector<Object*>::iterator object;
    … …  
    929929#undef __FUNCT__               
    930930#define __FUNCT__ "UpdateFromInputs"
    931 void  DataSet::UpdateFromInputs(ParameterInputs* inputs){
     931void  DataSet::UpdateFromInputs(void* inputs){
    932932
    933933        vector<Object*>::iterator object;
    … …  
    965965#undef __FUNCT__
    966966#define __FUNCT__ "DataSet::PenaltyCreateKMatrix"
    967 void  DataSet::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){
     967void  DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
    968968
    969969        vector<Object*>::iterator object;
    … …  
    983983#undef __FUNCT__
    984984#define __FUNCT__ "DataSet::PenaltyCreatePVector"
    985 void  DataSet::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){
     985void  DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
    986986
    987987        vector<Object*>::iterator object;
    … …  
    10091009#undef __FUNCT__
    10101010#define __FUNCT__ "RiftConstraints"
    1011 void  DataSet::RiftConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type){
     1011void  DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){
    10121012       
    10131013        throw ErrorException(__FUNCT__," not implemented yet!");
    … …  
    10261026#undef __FUNCT__
    10271027#define __FUNCT__ "MeltingConstraints"
    1028 void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type){
     1028void  DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){
    10291029
    10301030        throw ErrorException(__FUNCT__," not implemented yet!");
    … …  
    10651065
    10661066
    1067 void  DataSet::Du(Vec du_g,double*  u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){
     1067void  DataSet::Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type){
    10681068
    10691069
    … …  
    10831083}
    10841084
    1085 void  DataSet::Gradj(Vec grad_g,double*  u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type){
     1085void  DataSet::Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type){
    10861086
    10871087
    … …  
    11011101}               
    11021102               
    1103 void  DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){
     1103void  DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type){
    11041104
    11051105        double J=0;;

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

    r202 r246  
    1010#define DATASET_H_
    1111
    12 #include "../objects/Object.h"
    13 #include "../objects/ParameterInputs.h"
    1412#include <vector>
    1513#include "../toolkits/toolkits.h"
     14#include "../objects/objects.h"
    1615
    17 using namespace std;
    18 
    19 class DataSet {
     16class DataSet{
    2017       
    2118        private:
     19               
    2220                /*internals: */
    23                 vector<Object*> objects;
     21                std::vector<Object*> objects;
    2422               
    2523                /*type of dataset: */
    … …  
    6361                void  SetSorting(int* in_sorted_ids,int* in_id_offsets);
    6462                void  Sort();
    65                 void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs, int analysis_type);
    66                 void  CreatePVector(Vec pg,ParameterInputs* inputs, int analysis_type);
    67                 void  UpdateFromInputs(ParameterInputs* inputs);
    68                 void  PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type);
    69                 void  PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type);
     63                void  CreateKMatrix(Mat Kgg,void* inputs, int analysis_type);
     64                void  CreatePVector(Vec pg,void* inputs, int analysis_type);
     65                void  UpdateFromInputs(void* inputs);
     66                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
     67                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
    7068                int   RiftIsPresent();
    71                 void  RiftConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type);
     69                void  RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type);
    7270                int   MeltingIsPresent();
    73                 void  MeltingConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type);
     71                void  MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type);
    7472                DataSet* Copy(void);
    75                 void  Du(Vec du_g,double*  u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
    76                 void  Gradj(Vec grad_g,double*  u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type);
    77                 void  Misfit(double* pJ, double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
     73                void  Du(Vec du_g,double*  u_g,double* u_g_obs,void* inputs,int analysis_type);
     74                void  Gradj(Vec grad_g,double*  u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
     75                void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type);
    7876                void  VelocityExtrude(Vec ug,double* ug_serial);
    7977                int   DeleteObject(Object* object);
    8078
    81 
    82 
    8379};
    84 
    8580
    8681/*This routine cannot be object oriented, but need for demarshalling: */

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp

    r128 r246  
    5757int StokesEnum(void){                   return           4; }
    5858
     59/*inputs: */
     60int InputEnum(void){                    return         310; }
     61
     62
    5963
    6064/*functions on enums: */

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

    r128 r246  
    2121int LoadEnum(void);
    2222int RgbEnum(void);
     23int InputEnum(void);
    2324
    2425/*formulations: */

  • issm/trunk/src/c/Makefile.am

    r132 r246  
    4242                                        ./objects/Matpar.h\
    4343                                        ./objects/Matpar.cpp\
     44                                        ./objects/Input.h\
     45                                        ./objects/Input.cpp\
    4446                                        ./objects/ParameterInputs.h\
    4547                                        ./objects/ParameterInputs.cpp\
    … …  
    262264                                        ./objects/Matpar.h\
    263265                                        ./objects/Matpar.cpp\
     266                                        ./objects/Input.h\
     267                                        ./objects/Input.cpp\
    264268                                        ./objects/ParameterInputs.h\
    265269                                        ./objects/ParameterInputs.cpp\

  • issm/trunk/src/c/UpdateFromInputsx/UpdateFromInputsx.cpp

    r137 r246  
    2626        loads->Configure(elements, loads, nodes, materials);
    2727        nodes->Configure(elements, loads, nodes, materials);
    28        
     28
    2929        /*Update elements, nodes, loads and materials from inputs: */
    3030        elements->UpdateFromInputs(inputs);

  • issm/trunk/src/c/io/ParameterInputsInit.cpp

    r1 r246  
    99#endif
    1010
    11 #ifdef _SERIAL_
    1211
    1312#include "../objects/objects.h"
    … …  
    1817#undef __FUNCT__
    1918#define __FUNCT__ "ParameterInputsInit"
    20 void  ParameterInputsInit(ParameterInputs** pinputs, const mxArray* inputhandle){
    21                
    22         int dummy;
     19
     20#ifdef _SERIAL_
     21void ParameterInputs::Init( void* vinput_handle){
     22       
    2323        int i;
    2424        int nelems;
    25         int field_name_length;
     25        mxArray* input=NULL;
     26        const char* field_name=NULL;
     27        ConstDataHandle input_handle=NULL;
     28        int type;
     29
     30        char* string=NULL;
     31        char* name=NULL;
     32        double scalar;
     33        int    integer;
     34        int    ndof;
     35        int    numberofnodes;
     36        double* vector=NULL;
     37
     38        input_handle=(ConstDataHandle)vinput_handle;
    2639       
    27         /*output: */
    28         ParameterInputs* inputs=NULL;
     40        /*number of inputs: */
     41        nelems=mxGetNumberOfElements(input_handle);
    2942
    30         /*First allocate ParameterInputs: */
    31         inputs=NewParameterInputs();
     43        /*Go through ParameterInputs matlab structure, and recover each input object. Use it for fill
     44         * the ParameterInputs* "C" class object: */
    3245
    33         /*Now populate the structure: */
     46        for(i=0;i<nelems;i++){
    3447
    35         /*How many fields? :*/
    36         inputs->nfields = mxGetNumberOfFields(inputhandle); //tells us how many parameters we have
    37        
    38         /*Check that we don't have a multi-dimensional structure array: */
    39         nelems  = mxGetNumberOfElements(inputhandle);
    40        
    41         if(nelems!=1){
    42                 _printf_("%s%s\n",__FUNCT__," error message: input parameter structure must be of size 1!");
    43                 mexErrMsgTxt(" ");
    44         }
    45        
    46         /*Allocate fields and field names: */
    47         inputs->fields=(double**)xmalloc(inputs->nfields*sizeof(double*));
    48         inputs->field_names=(char**)xmalloc(inputs->nfields*sizeof(char*));
    49         inputs->fields_length=(int*)xmalloc(inputs->nfields*sizeof(int));
     48                /*Get i'th matlab input: */
     49                input=mxGetField(input_handle,i,"input");
    5050
    51         /*Initialize, we never know: */
    52         for(i=0;i<inputs->nfields;i++){
    53                 inputs->fields[i]=NULL;
    54                 inputs->field_names[i]=NULL;
    55                 inputs->fields_length[i]=0;
     51                /*This input has the following fields: name, type, string, integer, double,  and vector. Plus ndof and numberofnodes: */
     52                //first get type
     53                type=(int)*mxGetPr(mxGetField(input,0,"type"));
     54
     55                /*Get name: */
     56                name=mxArrayToString(mxGetField(input,0,"name"));
     57
     58                if (strcmp(name,"")==0)continue; //this input has no name! can't register it. keep going.
     59
     60                //string
     61                switch(type){
     62                        case 0:
     63                                //we have a string, recover it.
     64                                string=mxArrayToString(mxGetField(input,0,"string"));
     65                                //add to inputs.
     66                                this->Add(name,string);
     67                                xfree((void**)&string);
     68                                break;
     69                        case 1:
     70                                //integer
     71                                integer=(int)*(mxGetPr(mxGetField(input,0,"integer")));
     72                                this->Add(name,integer);
     73                                break;
     74                        case 2:
     75                                //double
     76                                scalar=(double)*(mxGetPr(mxGetField(input,0,"double")));
     77                                this->Add(name,scalar);
     78                                break;
     79                        case 3:
     80                                //double vector
     81                                vector=mxGetPr(mxGetField(input,0,"vector"));
     82                                ndof=(int)*(mxGetPr(mxGetField(input,0,"ndof")));
     83                                numberofnodes=(int)*(mxGetPr(mxGetField(input,0,"numberofnodes")));
     84                                this->Add(name,vector,ndof,numberofnodes);
     85                                break;
     86                        default:
     87                                throw ErrorException(__FUNCT__,exprintf("%s%i"," unknow parameter input type ",type));
     88                }
     89
    5690        }
    5791
    5892
    59         /*Go through input structure fields and populate our ParameterInputs object:*/
    60         for (i=0;i<inputs->nfields;i++){
    61                
    62                 /*Recover field value: */
    63                 FetchData((void**)&inputs->fields[i],&inputs->fields_length[i],&dummy,mxGetFieldByNumber( inputhandle,0,i),"Matrix","Mat");
    64 
    65                 /*Recover field name: */
    66                 const char* field_name=mxGetFieldNameByNumber(inputhandle,i);
    67                 field_name_length=strlen(field_name)+1;
    68                 inputs->field_names[i]=(char*)xmalloc(field_name_length*sizeof(char));
    69                 strcpy(inputs->field_names[i],field_name);
    70         }
    71                
    72         /*Assign output poiters:*/
    73         *pinputs=inputs;
    74 
    7593}
     94#else
     95void ParameterInputs::Init( void* vinput_handle){
     96}
    7697#endif //#ifdef _SERIAL-

  • issm/trunk/src/c/objects/Element.cpp

    r1 r246  
    55
    66
    7 #include "./Object.h"
    87#include "./Element.h"
    98#include "../EnumDefinitions/EnumDefinitions.h"
    10 #include "stdio.h"
    11 #include "../DataSet/DataSet.h"
    129
    1310int Element::Enum(void){

  • issm/trunk/src/c/objects/Element.h

    r1 r246  
    1010
    1111#include "./Object.h"
    12 #include "./ParameterInputs.h"
    1312#include "../toolkits/toolkits.h"
    14 #include "./Node.h"
    15 #include "./Matpar.h"
    1613
    1714class Element: public Object{
    … …  
    2724                virtual void  Demarshall(char** pmarshalled_dataset)=0;
    2825                virtual void  Configure(void* loads,void* nodes,void* materials)=0;
    29                 virtual void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type)=0;
    30                 virtual void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type)=0;
    31                 virtual void  UpdateFromInputs(ParameterInputs* inputs)=0;
    32                 virtual void  GetNodes(Node** nodes)=0;
    33                 virtual Matpar* GetMatPar()=0;
     26                virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0;
     27                virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type)=0;
     28                virtual void  UpdateFromInputs(void* inputs)=0;
     29                virtual void  GetNodes(void** nodes)=0;
     30                virtual void* GetMatPar()=0;
    3431                virtual int   GetShelf()=0;
    3532                virtual int   GetOnBed()=0;
    36                 virtual void          GetThicknessList(double* thickness_list)=0;
    37                 virtual void          GetBedList(double* bed_list)=0;
    38                 virtual void  Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type)=0;
    39                 virtual void  Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type)=0;
    40                 virtual void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type)=0;
    41                 virtual void  GradjB(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type)=0;
    42         virtual double Misfit(double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type)=0;
     33                virtual void  GetThicknessList(double* thickness_list)=0;
     34                virtual void  GetBedList(double* bed_list)=0;
     35                virtual void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0;
     36                virtual void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type)=0;
     37                virtual void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0;
     38                virtual void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0;
     39        virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0;
    4340
    4441                int           Enum();

  • issm/trunk/src/c/objects/FemModel.h

    r1 r246  
    99#include "../DataSet/DataSet.h"
    1010
    11 struct FemModel {
    12        
     11class DataSet;
     12
     13struct FemModel{
     14
    1315        DataSet*            elements;
    1416        DataSet*            nodes;

  • issm/trunk/src/c/objects/Icefront.cpp

    r181 r246  
    240240#define __FUNCT__ "Icefront::CreateKMatrix"
    241241
    242 void  Icefront::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     242void  Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
    243243
    244244        /*No stiffness loads applied, do nothing: */
    … …  
    250250#undef __FUNCT__
    251251#define __FUNCT__ "Icefront::CreatePVector"
    252 void  Icefront::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){
     252void  Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type){
    253253
    254254        /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
    … …  
    264264#undef __FUNCT__
    265265#define __FUNCT__ "Icefront::CreatePVectorDiagnosticHoriz"
    266 void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type){
     266void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type){
    267267
    268268        /*Branck on the type of icefront: */
    … …  
    278278#define __FUNCT__ "Icefront::CreatePVectorDiagnosticHorizSegment"
    279279
    280 void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,ParameterInputs* inputs, int analysis_type){
     280void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type){
    281281
    282282        int i,j;
    … …  
    306306        Matpar* matpar=NULL;
    307307        Node**  element_nodes=NULL;
     308        ParameterInputs* inputs=NULL;
     309
     310        inputs=(ParameterInputs*)vinputs;
    308311
    309312        /*Recover material and fill parameters: */
    310         matpar=element->GetMatPar();
     313        matpar=(Matpar*)element->GetMatPar();
    311314        fill=element->GetShelf();
    312315
    … …  
    323326        if(element_type==TriaEnum())element_nodes=(Node**)xmalloc(3*sizeof(Node*));
    324327        if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*));
    325         element->GetNodes(element_nodes);
     328        element->GetNodes((void**)element_nodes);
    326329
    327330        /*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */
    … …  
    389392#undef __FUNCT__
    390393#define __FUNCT__ "Icefont::CreatePVectorDiagnosticHorizQuad"
    391 void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,ParameterInputs* inputs, int analysis_type){
     394void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type){
    392395
    393396        int i,j;
    … …  
    413416        Matpar* matpar=NULL;
    414417        Node**  element_nodes=NULL;
     418        ParameterInputs* inputs=NULL;   
    415419
    416420        /*quad grids: */
    … …  
    427431        double v35[3];
    428432        double v45[3];
    429                
     433
     434        inputs=(ParameterInputs*)vinputs;
    430435
    431436        /*check icefront is associated to a pentaelem: */
    … …  
    434439        }
    435440        /*Recover material and fill parameters: */
    436         matpar=element->GetMatPar();
     441        matpar=(Matpar*)element->GetMatPar();
    437442        fill=element->GetShelf();
    438443
    … …  
    446451        //Identify which grids are comprised in the quad:
    447452        if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*));
    448         element->GetNodes(element_nodes);
     453        element->GetNodes((void**)element_nodes);
    449454
    450455        grid1=-1; grid2=-1; grid3=-1; grid4=-1;
    … …  
    549554#undef __FUNCT__
    550555#define __FUNCT__ "Icefront::UpdateFromInputs"
    551 void  Icefront::UpdateFromInputs(ParameterInputs* inputs){
     556void  Icefront::UpdateFromInputs(void* vinputs){
    552557
    553558        int i;
    … …  
    561566        double* h_param=NULL;
    562567        double* b_param=NULL;
    563 
    564         /*Get dof list: */
    565         nodes[0]->GetDofList(&doflist0[0],&numberofdofspernode);
    566         nodes[1]->GetDofList(&doflist1[0],&numberofdofspernode);
     568        int  dofs[1]={0};
     569
     570        ParameterInputs* inputs=NULL;   
     571
     572        inputs=(ParameterInputs*)vinputs;
     573
    567574        if(strcmp(type,"quad")==0){
    568                 nodes[2]->GetDofList(&doflist2[0],&numberofdofspernode);
    569                 nodes[3]->GetDofList(&doflist3[0],&numberofdofspernode);
    570         }
    571 
    572         /*Update internal data if inputs holds new values: */
    573         h_param=ParameterInputsRecover(inputs,"thickness");
    574         b_param=ParameterInputsRecover(inputs,"bed");
    575 
    576         if(h_param){
    577                 h[0]=h_param[doflist0[0]];
    578                 h[1]=h_param[doflist1[0]];
    579                 if(strcmp(type,"quad")==0){
    580                         h[2]=h_param[doflist2[0]];
    581                         h[3]=h_param[doflist3[0]];
    582                 }
    583         }
    584 
    585         if(b_param){
    586                 b[0]=b_param[doflist0[0]];
    587                 b[1]=b_param[doflist1[0]];
    588                 if(strcmp(type,"quad")==0){
    589                         b[2]=b_param[doflist2[0]];
    590                         b[3]=b_param[doflist3[0]];
    591                 }
    592         }
     575                inputs->Recover("thickness",&h[0],1,dofs,4,(void**)nodes);
     576                inputs->Recover("bed",&b[0],1,dofs,4,(void**)nodes);
     577        }
     578        else{
     579                inputs->Recover("thickness",&h[0],1,dofs,2,(void**)nodes);
     580                inputs->Recover("bed",&h[0],1,dofs,2,(void**)nodes);
     581        }
     582
    593583
    594584}
    … …  
    699689#undef __FUNCT__
    700690#define __FUNCT__ "Icefront::PenaltyCreateKMatrix"
    701 void  Icefront::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){
     691void  Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
    702692        /*do nothing: */
    703693}
    … …  
    705695#undef __FUNCT__
    706696#define __FUNCT__ "Icefront::PenaltyCreatePVector"
    707 void  Icefront::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){
     697void  Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
    708698        /*do nothing: */
    709699}

  • issm/trunk/src/c/objects/Icefront.h

    r94 r246  
    1414#define ICEFRONTSTRING 20 //max string length
    1515
     16class Element;
    1617class Icefront: public Load {
    1718
    … …  
    5657                void  DistributeNumDofs(int* numdofspernode,int analysis_type);
    5758                void  Configure(void* elements,void* nodes,void* materials);
    58                 void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    59                 void  UpdateFromInputs(ParameterInputs* inputs);
     59                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
     60                void  UpdateFromInputs(void* inputs);
    6061               
    61                 void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type);
    62                 void  CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type);
    63                 void  CreatePVectorDiagnosticHorizSegment( Vec pg,ParameterInputs* inputs, int analysis_type);
    64                 void  CreatePVectorDiagnosticHorizQuad( Vec pg,ParameterInputs* inputs, int analysis_type);
     62                void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     63                void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type);
     64                void  CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type);
     65                void  CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type);
    6566                void  GetDofList(int* doflist,int* pnumberofdofs);
    6667                void  SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill);
    6768                void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    6869                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill);
    69                 void  PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type);
    70                 void  PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type);
     70                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
     71                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
    7172                Object* copy();
    7273};

  • issm/trunk/src/c/objects/Load.cpp

    r1 r246  
    55
    66
    7 #include "./Object.h"
    87#include "./Load.h"
    98#include "../EnumDefinitions/EnumDefinitions.h"
    10 #include "stdio.h"
    11 #include "../DataSet/DataSet.h"
    129
    1310int Load::Enum(void){

  • issm/trunk/src/c/objects/Load.h

    r1 r246  
    1010
    1111#include "./Object.h"
    12 #include "./ParameterInputs.h"
    1312#include "../toolkits/toolkits.h"
    1413
    … …  
    2524                virtual void  Demarshall(char** pmarshalled_dataset)=0;
    2625                virtual void  Configure(void* elements,void* nodes,void* materials)=0;
    27                 virtual void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type)=0;
    28                 virtual void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type)=0;
    29                 virtual void  UpdateFromInputs(ParameterInputs* inputs)=0;
    30                 virtual void  PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type)=0;
    31                 virtual void  PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type)=0;
     26                virtual void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0;
     27                virtual void  CreatePVector(Vec pg, void* inputs, int analysis_type)=0;
     28                virtual void  UpdateFromInputs(void* inputs)=0;
     29                virtual void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type)=0;
     30                virtual void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type)=0;
    3231
    3332                int           Enum();

  • issm/trunk/src/c/objects/Material.cpp

    r1 r246  
    55
    66
    7 #include "./Object.h"
    87#include "./Material.h"
    98#include "../EnumDefinitions/EnumDefinitions.h"
    10 #include "stdio.h"
    11 #include "../DataSet/DataSet.h"
    129
    1310int Material::Enum(void){

  • issm/trunk/src/c/objects/Material.h

    r1 r246  
    88
    99#include "./Object.h"
    10 #include "./ParameterInputs.h"
    1110#include "../toolkits/toolkits.h"
    1211
    … …  
    2221                virtual char* GetName()=0;
    2322                virtual void  Demarshall(char** pmarshalled_dataset)=0;
    24                 virtual void  UpdateFromInputs(ParameterInputs* inputs)=0;
     23                virtual void  UpdateFromInputs(void* inputs)=0;
    2524                int           Enum();
    2625               

  • issm/trunk/src/c/objects/Matice.cpp

    r89 r246  
    107107#undef __FUNCT__
    108108#define __FUNCT__ "Matice::UpdateFromInputs"
    109 void  Matice::UpdateFromInputs(ParameterInputs* inputs){
     109void  Matice::UpdateFromInputs(void* inputs){
    110110       
    111111        //throw ErrorException(__FUNCT__," not supported yet!");

  • issm/trunk/src/c/objects/Matice.h

    r90 r246  
    3030                int   MyRank();
    3131                void  DistributeNumDofs(int* numdofspernode,int analysis_type);
    32                 void  UpdateFromInputs(ParameterInputs* inputs);
     32                void  UpdateFromInputs(void* inputs);
    3333                void  SetB(double B_param);
    3434                void  GetViscosity2d(double* pviscosity, double* pepsilon);

  • issm/trunk/src/c/objects/Matpar.cpp

    r111 r246  
    153153#undef __FUNCT__
    154154#define __FUNCT__ "Matpar::UpdateFromInputs"
    155 void  Matpar::UpdateFromInputs(ParameterInputs* inputs){
     155void  Matpar::UpdateFromInputs(void* inputs){
    156156       
    157157        //throw ErrorException(__FUNCT__," not supported yet!");

  • issm/trunk/src/c/objects/Matpar.h

    r111 r246  
    4040                int   MyRank();
    4141                void  DistributeNumDofs(int* numdofspernode,int analysis_type);
    42                 void  UpdateFromInputs(ParameterInputs* inputs);
     42                void  UpdateFromInputs(void* inputs);
    4343                double GetG();
    4444                double GetRhoIce();

  • issm/trunk/src/c/objects/Node.cpp

    r202 r246  
    437437#undef __FUNCT__
    438438#define __FUNCT__ "Node::UpdateFromInputs"
    439 void  Node::UpdateFromInputs(ParameterInputs* inputs){
    440        
    441         int doflist[MAXDOFSPERNODE*3];
    442         int numberofdofspernode;
    443 
    444         /*element: */
    445         double* x_param=NULL;
    446         double* y_param=NULL;
    447         double* z_param=NULL;
    448 
    449         /*Get dof list: */
    450         GetDofList(&doflist[0],&numberofdofspernode);
     439void  Node::UpdateFromInputs(void* vinputs){
     440       
     441        ParameterInputs* inputs=NULL;
     442        Node* node=this;
     443        int* dof={0};
     444
     445        /*Recover parameter inputs: */
     446        inputs=(ParameterInputs*)vinputs;
    451447
    452448        /*Update internal data if inputs holds new values: */
    453         x_param=ParameterInputsRecover(inputs,"x");
    454         y_param=ParameterInputsRecover(inputs,"y");
    455         z_param=ParameterInputsRecover(inputs,"z");
    456 
    457         if(x_param)x[0]=x_param[doflist[0]];
    458         if(y_param)x[1]=y_param[doflist[0]];
    459         if(y_param)x[2]=z_param[doflist[0]];
    460 
     449        inputs->Recover("x",&x[0],1,dof,1,(void**)&node);
     450        inputs->Recover("y",&x[1],1,dof,1,(void**)&node);
     451        inputs->Recover("z",&x[2],1,dof,1,(void**)&node);
     452
     453       
    461454}
    462455
    … …  
    494487
    495488}
     489               
     490int   Node::GetDoflist1(void){
     491        return doflist1[0];
     492}

  • issm/trunk/src/c/objects/Node.h

    r202 r246  
    88#include "./Object.h"
    99#include "../toolkits/toolkits.h"
    10 #include "./ParameterInputs.h"
    1110
    1211#define MAXDOFSPERNODE 4
    … …  
    6766                void  CreateVecSets(Vec pv_g,Vec pv_m,Vec pv_n,Vec pv_f,Vec pv_s);
    6867                void  GetDofList(int* outdoflist,int* pnumberofdofspernode);
     68                int   GetDoflist1(void);
    6969                double GetX();
    7070                double GetY();
    7171                double GetZ();
    7272                Object* copy();
    73                 void  UpdateFromInputs(ParameterInputs* inputs);
     73                void  UpdateFromInputs(void* inputs);
    7474                void  Configure(void* pnodes);
    7575                Node* GetUpperNode();

  • issm/trunk/src/c/objects/OptArgs.h

    r46 r246  
    2424#include "./FemModel.h"
    2525#include "./ParameterInputs.h"
     26class ParameterInputs;
    2627
    2728struct OptArgs{

  • issm/trunk/src/c/objects/ParameterInputs.cpp

    r1 r246  
    1 /*!\file:  ParameterInputs.cpp
    2  * \brief  routines for loading from matlab the parameter inputs to our non-linear iterations.
    3  */
    4 
    5 
    6 #include "./ParameterInputs.h"
     1/*!\file ParameterInputs.c
     2 * \brief: implementation of the ParameterInputs object
     3 */
     4
     5
     6#ifdef HAVE_CONFIG_H
     7        #include "config.h"
     8#else
     9#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     10#endif
     11
     12#include "stdio.h"
     13#include <string.h>
     14#include "../EnumDefinitions/EnumDefinitions.h"
    715#include "../shared/shared.h"
     16#include "../include/typedefs.h"
    817#include "../include/macros.h"
    9 
    10 /*--------------------------------------------------
    11         NewParameterInputs
    12   --------------------------------------------------*/
    13 
    14 ParameterInputs* NewParameterInputs(void) {
    15        
    16         ParameterInputs*        inputs=NULL;
    17 
    18         inputs=(ParameterInputs*)xmalloc(sizeof(ParameterInputs));
    19 
    20         inputs->nfields=0;
    21         inputs->fields=NULL;
    22         inputs->field_names=NULL;
    23         inputs->fields_length=NULL;
    24        
    25         return inputs;
    26 }
    27 
    28 
    29 /*--------------------------------------------------
    30         ParameterInputsEcho     
    31   --------------------------------------------------*/
    32 int ParameterInputsEcho(ParameterInputs* inputs,int longecho){
    33 
    34         int noerr=1;
    35         int i,j;
    36 
    37         /*short echo: */
    38         if(!longecho){
    39                 _printf_("field names, field pointers: \n");
    40                 for (i=0;i<inputs->nfields;i++){
    41                         _printf_("%s: %p\n",inputs->field_names[i],inputs->fields[i]);
    42                 }
    43         }
    44         /*long echo: */
    45         else{
    46                 for (i=0;i<inputs->nfields;i++){
    47                         double* values=inputs->fields[i];
    48                         _printf_("field names: %s\n",inputs->field_names[i]);
    49                         _printf_("field values: \n");
    50                         for(j=0;j<inputs->fields_length[i];j++){
    51                                 printf("  %g\n",values[j]);
    52                         }
    53                 }
    54         }
    55         return noerr;
    56 }
    57 
    58 /*--------------------------------------------------
    59         ParameterInputsRecover 
    60   --------------------------------------------------*/
    61 double* ParameterInputsRecover(ParameterInputs* inputs,char* field_name){
    62 
    63         double* param=NULL;
    64         int i;
    65        
    66         /*Find field_name input field, and return the value: */
    67         for(i=0;i<inputs->nfields;i++){
    68                 if (strcmp(inputs->field_names[i],field_name)==0){
    69                         param=inputs->fields[i];
    70                 }
    71         }
    72         return param;
    73 }
    74 
    75 /*--------------------------------------------------
    76         DeleteParameterInputs   
    77   --------------------------------------------------*/
    78 void DeleteParameterInputs(ParameterInputs** pinputs){
    79 
    80         ParameterInputs* inputs=NULL;
    81         int i;
    82         char*  field_name=NULL;
    83         double* field=NULL;
    84 
    85         #ifdef _PARALLEL_ //inputs operation is only done on the cluster, on the serial side, we leave the matlab memory manager handle inputs.
    86         inputs=*pinputs;
    87 
    88         if(inputs){
    89 
    90                 for(i=0;i<inputs->nfields;i++){
    91                         field_name=inputs->field_names[i];
    92                         xfree((void**)&field_name);
    93 
    94                         /*In parallel, erase: */
    95                         field=inputs->fields[i];
    96                         xfree((void**)&field);
    97                 }
    98 
    99                 xfree((void**)&inputs->field_names);
    100                 xfree((void**)&inputs->fields);
    101                 xfree((void**)&inputs->fields_length);
    102         }
    103         xfree((void**)pinputs);
    104         #else
    105         #endif
    106         return;
    107 }
    108 
    109 /*--------------------------------------------------
    110         ParameterInputsAddFromVec       
    111   --------------------------------------------------*/
    112 int  ParameterInputsAddFromVec(ParameterInputs* inputs,Vec vector,char* pfield_name){
    113        
    114         int     noerr=1;
    115         double* double_vector=NULL;
    116         double* double_vector2=NULL;
    117         int     vector_size;
    118         char*   field_name;
    119         char**  field_names;
    120         int     field_name_size;
    121         double** fields=NULL;
    122         int*    fields_length;
    123         int     i;
    124 
    125         /*Take vector, serialize it, and plug it into inputs with field name "field_name": */
    126         if (vector){
    127                 VecGetSize(vector,&vector_size);
    128                 VecToMPISerial(&double_vector,vector);
     18#include "./objects.h"
     19
     20ParameterInputs::ParameterInputs(){
     21
     22        dataset=new DataSet();
     23
     24}
    12925               
    130                 if(inputs->nfields==0){
    131                         /*Ok, we need to allocate fields, field_names and fields_length: */
    132                         inputs->nfields=1;
    133                        
    134                         inputs->field_names=(char**)xmalloc(sizeof(char*));
    135                         field_name_size=strlen(pfield_name)+1;
    136                         field_name=(char*)xmalloc(field_name_size*sizeof(char));
    137                         strcpy(field_name,pfield_name);
    138                         inputs->field_names[0]=field_name;
    139 
    140                         inputs->fields=(double**)xmalloc(sizeof(double*));
    141                         inputs->fields[0]=double_vector;
    142 
    143                         inputs->fields_length=(int*)xmalloc(sizeof(int));
    144                         inputs->fields_length[0]=vector_size;
    145                 }
    146                 else{
    147                         /*First, figure out if there is already an existing field with name pfield_name: */
    148                         for(i=0;i<inputs->nfields;i++){
    149                                 if (strcmp(inputs->field_names[i],pfield_name)==0){
    150                                         /*Erase old field and plug new field: */
    151                                         double_vector2=inputs->fields[i];
    152                                         xfree((void**)&double_vector2); //no leak!
    153                                         inputs->fields[i]=double_vector;
    154                                         inputs->fields_length[i]=vector_size;
    155                                         return noerr;
    156                                 }
    157                         }
    158 
    159                         /*Ok, if we are here, the field did not already exist, we need to add it. First reallocate
    160                          * the inputs: */
    161                         field_names=inputs->field_names; //save field_names pointer.
    162                         inputs->field_names=(char**)xmalloc((inputs->nfields+1)*sizeof(char*));
    163                         for(i=0;i<inputs->nfields;i++){
    164                                 inputs->field_names[i]=field_names[i];
    165                         }
    166                         xfree((void**)&field_names);
    167 
    168                         field_name_size=strlen(pfield_name)+1;
    169                         field_name=(char*)xmalloc(field_name_size*sizeof(char));
    170                         strcpy(field_name,pfield_name);
    171                         inputs->field_names[inputs->nfields]=field_name;
    172 
    173 
    174 
    175                         fields=inputs->fields; //save fields pointer
    176                         inputs->fields=(double**)xmalloc((inputs->nfields+1)*sizeof(double*));
    177                         for(i=0;i<inputs->nfields;i++){
    178                                 inputs->fields[i]=fields[i];
    179                         }
    180                         xfree((void**)&fields);
    181                         inputs->fields[inputs->nfields]=double_vector;
    182 
    183                         fields_length=inputs->fields_length;
    184                         inputs->fields_length=(int*)xmalloc((inputs->nfields+1)*sizeof(int));
    185                         for(i=0;i<inputs->nfields;i++){
    186                                 inputs->fields_length[i]=fields_length[i];
    187                         }
    188                         inputs->fields_length[inputs->nfields]=vector_size;
    189                         xfree((void**)&fields_length);
    190                        
    191                         inputs->nfields++;
    192                 }
    193         }
    194         else{
    195                 /*Do nothing, there is no input to add!: */
    196                 ;
    197         }
    198         return noerr;
    199 }
    200 
    201 /*--------------------------------------------------
    202         ParameterInputsAddFromDouble   
    203   --------------------------------------------------*/
    204 int  ParameterInputsAddFromDouble(ParameterInputs* inputs,double scalar, char* pfield_name){
    205        
    206         int     noerr=1;
    207         double* double_vector2=NULL;
    208         double* double_vector=NULL;
    209         int     vector_size;
    210         char*   field_name;
    211         char**  field_names;
    212         int     field_name_size;
    213         double** fields=NULL;
    214         int*    fields_length;
    215         int     i;
    216 
    217         /*Create double_vector pointing to scalar: */
    218         vector_size=1;
    219         double_vector=(double*)xmalloc(sizeof(double));
    220         double_vector[0]=scalar;
    221 
    222 
    223         if(inputs->nfields==0){
    224                 /*Ok, we need to allocate fields, field_names and fields_length: */
    225                 inputs->nfields=1;
     26ParameterInputs::~ParameterInputs(){
     27
     28        delete dataset;
     29
     30}
     31
     32void ParameterInputs::Echo(){
     33
     34        printf("ParameterInputs echo: \n");
     35        dataset->Echo();
     36
     37}
     38
     39void ParameterInputs::purge(char* name){
     40
     41        int i;
     42        Input* input=NULL;
     43
     44        /*Go through dataset, and figure out if an Input
     45         * already has the name "name". If so, delete it : */
     46       
     47        for(i=0;i<dataset->Size();i++){
     48                input=(Input*)dataset->GetObjectByOffset(i);
     49
     50                if (input->IsSameName(name)){
     51                        /*delete object: */
     52                        dataset->DeleteObject(input);
     53                }
     54        }
     55}
    22656               
    227                 inputs->field_names=(char**)xmalloc(sizeof(char*));
    228                 field_name_size=strlen(pfield_name)+1;
    229                 field_name=(char*)xmalloc(field_name_size*sizeof(char));
    230                 strcpy(field_name,pfield_name);
    231                 inputs->field_names[0]=field_name;
    232 
    233                 inputs->fields=(double**)xmalloc(sizeof(double*));
    234                 inputs->fields[0]=double_vector;
    235 
    236                 inputs->fields_length=(int*)xmalloc(sizeof(int));
    237                 inputs->fields_length[0]=vector_size;
    238         }
    239         else{
    240                 /*First, figure out if there is already an existing field with name pfield_name: */
    241                 for(i=0;i<inputs->nfields;i++){
    242                         if (strcmp(inputs->field_names[i],pfield_name)==0){
    243                                 /*Erase old field and plug new field: */
    244                                 double_vector2=inputs->fields[i];
    245                                 xfree((void**)&double_vector2); //no leak!
    246                                 inputs->fields[i]=double_vector;
    247                                 inputs->fields_length[i]=vector_size;
    248                                 return noerr;
    249                         }
    250                 }
    251 
    252                 /*Ok, if we are here, the field did not already exist, we need to add it. First reallocate
    253                  * the inputs: */
    254                 field_names=inputs->field_names; //save field_names pointer.
    255                 inputs->field_names=(char**)xmalloc((inputs->nfields+1)*sizeof(char*));
    256                 for(i=0;i<inputs->nfields;i++){
    257                         inputs->field_names[i]=field_names[i];
    258                 }
    259                 xfree((void**)&field_names);
    260 
    261                 field_name_size=strlen(pfield_name)+1;
    262                 field_name=(char*)xmalloc(field_name_size*sizeof(char));
    263                 strcpy(field_name,pfield_name);
    264                 inputs->field_names[inputs->nfields]=field_name;
    265 
    266 
    267 
    268                 fields=inputs->fields; //save fields pointer
    269                 inputs->fields=(double**)xmalloc((inputs->nfields+1)*sizeof(double*));
    270                 for(i=0;i<inputs->nfields;i++){
    271                         inputs->fields[i]=fields[i];
    272                 }
    273                 xfree((void**)&fields);
    274                 inputs->fields[inputs->nfields]=double_vector;
    275 
    276                 fields_length=inputs->fields_length;
    277                 inputs->fields_length=(int*)xmalloc((inputs->nfields+1)*sizeof(int));
    278                 for(i=0;i<inputs->nfields;i++){
    279                         inputs->fields_length[i]=fields_length[i];
    280                 }
    281                 inputs->fields_length[inputs->nfields]=vector_size;
    282                 xfree((void**)&fields_length);
    283                
    284                 inputs->nfields++;
    285         }
    286         return noerr;
    287 }
    288 
    289 /*--------------------------------------------------
    290         ParameterInputsAddFromMat       
    291   --------------------------------------------------*/
    292 int  ParameterInputsAddFromMat(ParameterInputs* inputs,double* vector,int vector_size,char* pfield_name){
    293        
    294         int     noerr=1;
    295         char*   field_name;
    296         char**  field_names;
    297         int     field_name_size;
    298         double** fields=NULL;
    299         int*    fields_length;
    300         double* double_vector=NULL;
    301         double* double_vector2=NULL;
    302         int     i;
    303 
    304         /*Make a copy of input vector, so as to not create any leaks: */
    305         double_vector=(double*)xmalloc(vector_size*sizeof(double));
    306         memcpy(double_vector,vector,vector_size*sizeof(double));
    307 
    308         /*Take vector, serialize it, and plug it into inputs with field name "field_name": */
    309         if (vector){
    310                
    311                 if(inputs->nfields==0){
    312                         /*Ok, we need to allocate fields, field_names and fields_length: */
    313                         inputs->nfields=1;
    314                        
    315                         inputs->field_names=(char**)xmalloc(sizeof(char*));
    316                         field_name_size=strlen(pfield_name)+1;
    317                         field_name=(char*)xmalloc(field_name_size*sizeof(char));
    318                         strcpy(field_name,pfield_name);
    319                         inputs->field_names[0]=field_name;
    320 
    321                         inputs->fields=(double**)xmalloc(sizeof(double*));
    322                         inputs->fields[0]=double_vector;
    323 
    324                         inputs->fields_length=(int*)xmalloc(sizeof(int));
    325                         inputs->fields_length[0]=vector_size;
    326                 }
    327                 else{
    328 
    329                         /*First, figure out if there is already an existing field with name pfield_name: */
    330                         for(i=0;i<inputs->nfields;i++){
    331                                 if (strcmp(inputs->field_names[i],pfield_name)==0){
    332                                         /*Just xfree old field and plug new field: */
    333                                         double_vector2=inputs->fields[i];
    334                                         xfree((void**)&double_vector2);
    335                                         inputs->fields[i]=double_vector;
    336                                         inputs->fields_length[i]=vector_size;
    337                                         return noerr;
    338                                 }
    339                         }
    340 
    341                         /*Ok, if we are here, the field did not already exist, we need to add it. First reallocate
    342                          * the inputs: */
    343                         field_names=inputs->field_names; //save field_names pointer.
    344                         inputs->field_names=(char**)xmalloc((inputs->nfields+1)*sizeof(char*));
    345                         for(i=0;i<inputs->nfields;i++){
    346                                 inputs->field_names[i]=field_names[i];
    347                         }
    348                         xfree((void**)&field_names);
    349 
    350                         field_name_size=strlen(pfield_name)+1;
    351                         field_name=(char*)xmalloc(field_name_size*sizeof(char));
    352                         strcpy(field_name,pfield_name);
    353                         inputs->field_names[inputs->nfields]=field_name;
    354 
    355 
    356                         fields=inputs->fields; //save fields pointer
    357                         inputs->fields=(double**)xmalloc((inputs->nfields+1)*sizeof(double*));
    358                         for(i=0;i<inputs->nfields;i++){
    359                                 inputs->fields[i]=fields[i];
    360                         }
    361                         xfree((void**)&fields);
    362                         inputs->fields[inputs->nfields]=double_vector;
    363 
    364                         fields_length=inputs->fields_length;
    365                         inputs->fields_length=(int*)xmalloc((inputs->nfields+1)*sizeof(int));
    366                         for(i=0;i<inputs->nfields;i++){
    367                                 inputs->fields_length[i]=fields_length[i];
    368                         }
    369                         inputs->fields_length[inputs->nfields]=vector_size;
    370                         xfree((void**)&fields_length);
    371                        
    372                         inputs->nfields++;
    373                 }
    374         }
    375         else{
    376                 /*Do nothing, there is no input to add!: */
    377                 ;
    378         }
    379         return noerr;
    380 }
     57
     58void ParameterInputs::Add(char* name,char* string){
     59
     60        Input* input=NULL;
     61
     62        /*First, purge object with same name: */
     63        this->purge(name);
     64
     65        /*We are sure an input of the same name does not exist. Create new
     66         * input: */
     67        input=new Input(name,string);
     68
     69        /*Add input to dataset: */
     70        dataset->AddObject(input);
     71
     72}
     73
     74void ParameterInputs::Add(char* name,int integer){
     75       
     76        Input* input=NULL;
     77
     78        /*First, purge object with same name: */
     79        this->purge(name);
     80
     81        /*We are sure an input of the same name does not exist. Create new
     82         * input: */
     83        input=new Input(name,integer);
     84
     85        /*Add input to dataset: */
     86        dataset->AddObject(input);
     87}
     88
     89
     90void ParameterInputs::Add(char* name,double scalar){
     91
     92        Input* input=NULL;
     93
     94        /*First, purge object with same name: */
     95        this->purge(name);
     96
     97        /*We are sure an input of the same name does not exist. Create new
     98         * input: */
     99        input=new Input(name,scalar);
     100
     101        /*Add input to dataset: */
     102        dataset->AddObject(input);
     103
     104}
     105
     106void ParameterInputs::Add(char* name,double* vector,int ndof,int numberofnodes){
     107       
     108        Input* input=NULL;
     109
     110        /*First, purge object with same name: */
     111        this->purge(name);
     112
     113        /*We are sure an input of the same name does not exist. Create new
     114         * input: */
     115        input=new Input(name,vector,ndof,numberofnodes);
     116
     117        /*Add input to dataset: */
     118        dataset->AddObject(input);
     119
     120}
     121
     122void ParameterInputs::Add(char* name,Vec petscvector,int ndof, int numberofnodes){
     123       
     124        Input* input=NULL;
     125
     126        /*First, purge object with same name: */
     127        this->purge(name);
     128
     129        /*We are sure an input of the same name does not exist. Create new
     130         * input: */
     131        input=new Input(name,petscvector,ndof,numberofnodes);
     132
     133        /*Add input to dataset: */
     134        dataset->AddObject(input);
     135
     136}
     137int ParameterInputs::Recover(char* name, char** pstring){
     138
     139        int found=0;
     140        int i;
     141        Input* input=NULL;
     142        char* string=NULL;
     143
     144        /*Go through dataset, and figure out if an Input
     145         * has the name "name": */
     146        for(i=0;i<dataset->Size();i++){
     147                input=(Input*)dataset->GetObjectByOffset(i);
     148
     149                if (input->IsSameName(name)){
     150
     151                        /*Get string out of this input: */
     152                        input->Recover(&string);
     153                        found=1;break;
     154                }
     155        }
     156        /*Assign output pointers:*/
     157        *pstring=string;
     158        return found;
     159
     160}
     161int ParameterInputs::Recover(char* name, int* pinteger){
     162       
     163        int found=0;
     164        int i;
     165        Input* input=NULL;
     166        int integer;
     167
     168        /*Go through dataset, and figure out if an Input
     169         * has the name "name": */
     170        for(i=0;i<dataset->Size();i++){
     171                input=(Input*)dataset->GetObjectByOffset(i);
     172
     173                if (input->IsSameName(name)){
     174
     175                        /*Get string out of this input: */
     176                        input->Recover(&integer);
     177                        found=1; break;
     178
     179                }
     180        }
     181        /*Assign output pointers:*/
     182        *pinteger=integer;
     183
     184        return found;
     185}
     186int ParameterInputs::Recover(char* name, double* pscalar){
     187
     188        int found=0;
     189        int i;
     190        Input* input=NULL;
     191        double scalar;
     192
     193        /*Go through dataset, and figure out if an Input
     194         * has the name "name": */
     195        for(i=0;i<dataset->Size();i++){
     196                input=(Input*)dataset->GetObjectByOffset(i);
     197
     198                if (input->IsSameName(name)){
     199
     200                        /*Get string out of this input: */
     201                        input->Recover(&scalar);
     202                        found=1; break;
     203
     204                }
     205        }
     206        /*Assign output pointers:*/
     207        *pscalar=scalar;
     208
     209        return found;
     210
     211}
     212
     213#undef __FUNCT__
     214#define __FUNCT__ "ParameterInputs::Recover"
     215int ParameterInputs::Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** pnodes){
     216
     217        int i;
     218        Node** nodes=NULL;
     219        Input* input=NULL;
     220        int found=0;
     221
     222        /*First recover nodes pointer :*/
     223        nodes=(Node**)pnodes;
     224
     225        /*Go through dataset, and figure out if an Input
     226         * has the name "name": */
     227        for(i=0;i<dataset->Size();i++){
     228                input=(Input*)dataset->GetObjectByOffset(i);
     229
     230                if (input->IsSameName(name)){
     231                        found=1;
     232                        break;
     233                }
     234        }
     235
     236        if(found==0)return 0;
     237       
     238        /*call submethod: */
     239        input->Recover(values,ndof,dofs,numnodes,(void**)nodes);
     240        return 1;
     241}

  • issm/trunk/src/c/objects/ParameterInputs.h

    r1 r246  
    11/*!\file:  ParameterInputs.h
    2  * \brief prototype for parameter inputs, from matlab workspace.
     2 * \brief dataset filled with input objects. Used by CreateKMatrix and CreatePVector  to hold
     3 * parameters coming from the workspace.
    34 */
    4 #ifndef PARAMETER_INPUTS_H_
    5 #define PARAMETER_INPUTS_H_
     5
     6#ifndef _PARAMETERINPUTS_H_
     7#define _PARAMETERINPUTS_H_
    68
    79#include "../toolkits/toolkits.h"
     10#include "../DataSet/DataSet.h"
    811
    9 struct ParameterInputs{
     12class DataSet;
    1013
    11         int     nfields; //number of inputs
    12         char** field_names;  //characterization of each field.
    13         double** fields;  //value of each field.
    14         int* fields_length; //length of each values field
     14class ParameterInputs{
    1515
     16        private:
     17                DataSet* dataset; //cannot declare DataSet* dataset, to avoid circular referencing.
     18
     19        public:
     20                ParameterInputs();
     21                ~ParameterInputs();
     22                void Echo();
     23               
     24                void Add(char* name,char* string);
     25                void Add(char* name,int integer);
     26                void Add(char* name,double scalar);
     27                void Add(char* name,double* vector,int ndof,int numberofnodes);
     28                void Add(char* name,Vec petscvector,int ndof, int numberofnodes);
     29
     30                int  Recover(char* name, char** pstring);
     31                int  Recover(char* name, int* pinteger);
     32                int  Recover(char* name, double* pscalar);
     33                int  Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** nodes); //void** because otherwise we get a circular dependency with Node
     34
     35                /*declaration found in io: */
     36                void Init( void* data_handle);
     37                void purge(char* name);
    1638};
    1739
    18 ParameterInputs* NewParameterInputs();
    19 void             DeleteParameterInputs(ParameterInputs** pinputs);
    20 int              ParameterInputsEcho(ParameterInputs* inputs,int longecho);
    21 double*          ParameterInputsRecover(ParameterInputs* inputs,char* field_name);
    22 
    23 int              ParameterInputsAddFromVec(ParameterInputs* inputs,Vec vector,char* field_name);
    24 int              ParameterInputsAddFromMat(ParameterInputs* inputs,double* vector,int vector_size,char* field_name);
    25 int              ParameterInputsAddFromDouble(ParameterInputs* inputs,double scalar, char* pfield_name);
    26 
    27 #ifdef _SERIAL_
    28 #include "mex.h"
    29 void  ParameterInputsInit(ParameterInputs** pthis, const mxArray* inputhandle);
    3040#endif
    31 
    32 #endif

  • issm/trunk/src/c/objects/Penpair.cpp

    r234 r246  
    208208#define __FUNCT__ "Penpair::CreateKMatrix"
    209209
    210 void  Penpair::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     210void  Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
    211211
    212212        /*No loads applied, do nothing: */
    … …  
    217217#undef __FUNCT__
    218218#define __FUNCT__ "Penpair::CreatePVector"
    219 void  Penpair::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){
     219void  Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type){
    220220
    221221        /*No loads applied, do nothing: */
    … …  
    225225#undef __FUNCT__
    226226#define __FUNCT__ "Penpair::UpdateFromInputs"
    227 void  Penpair::UpdateFromInputs(ParameterInputs* inputs){
     227void  Penpair::UpdateFromInputs(void* inputs){
    228228       
    229229}
    … …  
    231231#undef __FUNCT__
    232232#define __FUNCT__ "Penpair::PenaltyCreateKMatrix"
    233 void  Penpair::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){
     233void  Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
    234234        if(numdofs==1){
    235235               
    … …  
    245245#undef __FUNCT__
    246246#define __FUNCT__ "Penpair::PenaltyCreatePVector"
    247 void  Penpair::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){
     247void  Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
    248248        if(numdofs==1){
    249249               

  • issm/trunk/src/c/objects/Penpair.h

    r1 r246  
    99#include "./Element.h"
    1010
     11class Element;
    1112class Penpair: public Load{
    1213
    … …  
    5253                void  DistributeNumDofs(int* numdofspernode,int analysis_type);
    5354                void  Configure(void* elements,void* nodes,void* materials);
    54                 void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    55                 void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type);
    56                 void  UpdateFromInputs(ParameterInputs* inputs);
    57                 void  PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type);
    58                 void  PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type);
     55                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
     56                void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     57                void  UpdateFromInputs(void* inputs);
     58                void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
     59                void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
    5960                Object* copy();
    6061};

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

    r241 r246  
    279279#define __FUNCT__ "Penta::CreateKMatrix"
    280280
    281 void  Penta::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     281void  Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
    282282
    283283        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    … …  
    301301#undef __FUNCT__
    302302#define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz"
    303 void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, ParameterInputs* inputs, int analysis_type){
     303void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type){
    304304
    305305
    … …  
    368368
    369369        /*input parameters for structural analysis (diagnostic): */
    370         double* velocity_param=NULL;
    371         double  vxvy_list[numgrids][2];
    372         double* oldvelocity_param=NULL;
    373         double  oldvxvy_list[numgrids][2];
     370        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     371        double  oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     372        int     dofs[2]={0,1};
    374373        double  thickness;
    375374
    … …  
    381380        double MOUNTAINKEXPONENT=10;
    382381
     382        ParameterInputs* inputs=NULL;
     383
    383384        /*Collapsed formulation: */
    384385        Tria*  tria=NULL;
     386
     387        /*recover pointers: */
     388        inputs=(ParameterInputs*)vinputs;
    385389
    386390        /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the
    … …  
    397401                /*This element should be collapsed into a tria element at its base. Create this tria element,
    398402                 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
    399                 tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     403                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    400404                tria->CreateKMatrix(Kgg,inputs, analysis_type);
    401405                delete tria;
    … …  
    407411
    408412                /*recover extra inputs from users, at current convergence iteration: */
    409                 velocity_param=ParameterInputsRecover(inputs,"velocity");
    410                 oldvelocity_param=ParameterInputsRecover(inputs,"old_velocity");
     413                inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     414                inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    411415
    412416                /* Get node coordinates and dof list: */
    … …  
    418422                        for(j=0;j<numdofs;j++){
    419423                                Ke_gg[i][j]=0.0;
    420                         }
    421                 }
    422 
    423                 /*Initialize vxvy_list: */
    424                 for(i=0;i<numgrids;i++){
    425                         if(velocity_param){
    426                                 vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
    427                                 vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
    428                         }
    429                         else{
    430                                 vxvy_list[i][0]=0;
    431                                 vxvy_list[i][1]=0;
    432                         }
    433                 }
    434                        
    435                 /*Initialize oldvxvy_list: */
    436                 for(i=0;i<numgrids;i++){
    437                         if(oldvelocity_param){
    438                                 oldvxvy_list[i][0]=oldvelocity_param[doflist[i*numberofdofspernode+0]];
    439                                 oldvxvy_list[i][1]=oldvelocity_param[doflist[i*numberofdofspernode+1]];
    440                         }
    441                         else{
    442                                 oldvxvy_list[i][0]=0;
    443                                 oldvxvy_list[i][1]=0;
    444424                        }
    445425                }
    … …  
    552532                         * grids: */
    553533
    554                         tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     534                        tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    555535                        tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type);
    556536                        delete tria;
    … …  
    574554#undef __FUNCT__
    575555#define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert"
    576 void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type){
     556void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type){
    577557
    578558        /* local declarations */
    … …  
    610590        double  Bprime[NDOF1][numgrids];
    611591        double  DL_scalar;
    612        
     592
     593        ParameterInputs* inputs=NULL;
     594
    613595        /*Collapsed formulation: */
    614596        Tria*  tria=NULL;
     597
     598        /*recover pointers: */
     599        inputs=(ParameterInputs*)vinputs;
     600       
    615601
    616602        /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness
    617603         * matrix: */
    618604        if(onsurface){
    619                 tria=SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
     605                tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface
    620606                tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type);
    621607                delete tria;
    … …  
    706692#undef __FUNCT__
    707693#define __FUNCT__ "Penta::CreatePVector"
    708 void  Penta::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){
     694void  Penta::CreatePVector(Vec pg, void* inputs, int analysis_type){
    709695
    710696        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    … …  
    724710#undef __FUNCT__
    725711#define __FUNCT__ "Penta::UpdateFromInputs"
    726 void  Penta::UpdateFromInputs(ParameterInputs* inputs){
    727 
    728         int i;
    729         int doflist[MAXDOFSPERNODE*3];
    730         int numberofdofspernode;
    731 
    732         /*element: */
    733         double* h_param=NULL;
    734         double* s_param=NULL;
    735         double* b_param=NULL;
    736         double* k_param=NULL;
    737         double* melting_param=NULL;
    738         double* accumulation_param=NULL;
    739        
    740         /*material: */
    741         double* temperature_average_param=NULL;
     712void  Penta::UpdateFromInputs(void* vinputs){
     713
     714        int     dofs[1]={0};
     715        double  temperature_list[6];
    742716        double  temperature_average;
    743         double* B_param=NULL;
     717        double  B_list[6];
    744718        double  B_average;
    745719
    746 
    747         /*Get dof list: */
    748         GetDofList(&doflist[0],&numberofdofspernode);
     720        ParameterInputs* inputs=NULL;
     721
     722        /*recover pointers: */
     723        inputs=(ParameterInputs*)vinputs;
    749724
    750725        /*Update internal data if inputs holds new values: */
    751         h_param=ParameterInputsRecover(inputs,"thickness");
    752         s_param=ParameterInputsRecover(inputs,"surface");
    753         b_param=ParameterInputsRecover(inputs,"bed");
    754         k_param=ParameterInputsRecover(inputs,"drag");
    755         melting_param=ParameterInputsRecover(inputs,"melting");
    756         accumulation_param=ParameterInputsRecover(inputs,"accumulation");
    757 
    758         for(i=0;i<6;i++){
    759                 if(h_param)h[i]=h_param[doflist[i*numberofdofspernode+0]];
    760                 if(s_param)s[i]=s_param[doflist[i*numberofdofspernode+0]];
    761                 if(b_param)b[i]=b_param[doflist[i*numberofdofspernode+0]];
    762                 if(k_param)k[i]=k_param[doflist[i*numberofdofspernode+0]];
    763                 if(melting_param)melting[i]=melting_param[doflist[i*numberofdofspernode+0]];
    764                 if(accumulation_param)accumulation[i]=accumulation_param[doflist[i*numberofdofspernode+0]];
    765         }
    766 
     726        inputs->Recover("thickness",&h[0],1,dofs,6,(void**)nodes);
     727        inputs->Recover("surface",&s[0],1,dofs,6,(void**)nodes);
     728        inputs->Recover("bed",&b[0],1,dofs,6,(void**)nodes);
     729        inputs->Recover("drag",&k[0],1,dofs,6,(void**)nodes);
     730        inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes);
     731        inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes);
     732       
    767733        //Update material if necessary
    768         temperature_average_param=ParameterInputsRecover(inputs,"temperature_average");
    769         B_param=ParameterInputsRecover(inputs,"B");
    770         if(temperature_average_param){
    771                 temperature_average=0;
    772                 for(i=0;i<6;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
    773                 temperature_average=temperature_average/6.0;
     734        if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){
     735                temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0;
    774736                B_average=Paterson(temperature_average);
    775737                matice->SetB(B_average);
    776738        }
    777 
    778         if(B_param){
    779                 B_average=0;
    780                 for(i=0;i<6;i++) B_average+=B_param[doflist[i*numberofdofspernode+0]];
    781                 B_average=B_average/6.0;
     739       
     740        if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){
     741                B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0;
    782742                matice->SetB(B_average);
    783743        }
    … …  
    785745}
    786746
    787 Matpar* Penta::GetMatPar(){
     747void* Penta::GetMatPar(){
    788748        return matpar;
    789749}
    … …  
    794754
    795755
    796 void  Penta::GetNodes(Node** pnodes){
     756void  Penta::GetNodes(void** vpnodes){
    797757        int i;
     758        Node** pnodes=(Node**)vpnodes;
     759
    798760        for(i=0;i<6;i++){
    799761                pnodes[i]=nodes[i];
    … …  
    823785#undef __FUNCT__
    824786#define __FUNCT__ "Penta::Du"
    825 void  Penta::Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){
     787void  Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type){
    826788
    827789        int i;
    … …  
    834796        else{
    835797               
    836                 tria=SpawnTria(3,4,5); //grids 0, 1 and 2 make the new tria.
     798                tria=(Tria*)SpawnTria(3,4,5); //grids 0, 1 and 2 make the new tria.
    837799                tria->Du(du_g,u_g,u_g_obs,inputs,analysis_type);
    838800                delete tria;
    … …  
    843805#undef __FUNCT__
    844806#define __FUNCT__ "Penta::Gradj"
    845 void  Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type,char* control_type){
     807void  Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){
    846808
    847809        if (strcmp(control_type,"drag")==0){
    … …  
    856818#undef __FUNCT__
    857819#define __FUNCT__ "Penta::GradjDrag"
    858 void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type){
     820void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){
    859821       
    860822       
    … …  
    867829        else{
    868830               
    869                 tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     831                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    870832                tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type);
    871833                delete tria;
    … …  
    876838#undef __FUNCT__
    877839#define __FUNCT__ "Penta::GradjB"
    878 void  Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type){
     840void  Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){
    879841        throw ErrorException(__FUNCT__," not supported yet!");
    880842}
    … …  
    882844#undef __FUNCT__
    883845#define __FUNCT__ "Penta::Misfit"
    884 double Penta::Misfit(double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){
     846double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type){
    885847       
    886848        double J;
    … …  
    894856
    895857                /*use the Tria::CreateMisfit routine, on a Tria Element which is made of the 3 upper grids at the surface of the ice sheet. : */
    896                 tria=SpawnTria(3,4,5);
     858                tria=(Tria*)SpawnTria(3,4,5);
    897859                J=tria->Misfit( velocity,obs_velocity,inputs,analysis_type);
    898860                delete tria;
    … …  
    903865#undef __FUNCT__
    904866#define __FUNCT__ "Penta::SpawnTria"
    905 Tria*  Penta::SpawnTria(int g0, int g1, int g2){
     867void*  Penta::SpawnTria(int g0, int g1, int g2){
    906868
    907869        /*out of grids g0,g1 and g2 from Penta, build a tria element: */
    … …  
    13451307#undef __FUNCT__
    13461308#define __FUNCT__ "Penta::CreatePVectorDiagnosticHoriz"
    1347 void Penta::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs,int analysis_type){
     1309void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type){
    13481310
    13491311        int i,j;
    … …  
    13881350        double  pe_g_gaussian[numdofs];
    13891351
     1352        ParameterInputs* inputs=NULL;
     1353
    13901354        /*Spawning: */
    13911355        Tria* tria=NULL;
     1356
     1357        /*recover pointers: */
     1358        inputs=(ParameterInputs*)vinputs;
    13921359
    13931360
    … …  
    14051372                /*This element should be collapsed into a tria element at its base. Create this tria element,
    14061373                 *and use its CreatePVector functionality to return an elementary load vector: */
    1407                 tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     1374                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    14081375                tria->CreatePVector(pg,inputs, analysis_type);
    14091376                delete tria;
    … …  
    16441611#undef __FUNCT__
    16451612#define __FUNCT__ "Penta:CreatePVectorDiagnosticVert"
    1646 void  Penta::CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type){
     1613void  Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type){
    16471614
    16481615        int i;
    … …  
    16851652
    16861653        /*input parameters for structural analysis (diagnostic): */
    1687         double* velocity_param=NULL;
    1688         double vx_list[numgrids];
    1689         double vy_list[numgrids];
     1654        double vx_list[numgrids]={0,0,0};
     1655        double vy_list[numgrids]={0,0,0};
    16901656        double du[3];
    16911657        double dv[3];
    16921658        double dudx,dvdy;
     1659        int     dofs1[1]={0};
     1660        int     dofs2[1]={1};
     1661
     1662        ParameterInputs* inputs=NULL;
     1663
     1664        /*recover pointers: */
     1665        inputs=(ParameterInputs*)vinputs;
    16931666
    16941667        /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the
    16951668         *diagnostic base vertical stifness: */
    16961669        if(onbed){
    1697                 tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
     1670                tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock
    16981671                tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type);
    16991672                delete tria;
    … …  
    17091682
    17101683        /* recover input parameters: */
    1711         velocity_param=ParameterInputsRecover(inputs,"velocity");
    1712 
    1713         if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
    1714 
    1715         /*Initialize vx_list and vy_list: */
    1716         for(i=0;i<numgrids;i++){
    1717                 /*doflist is a list of dofs for the node. We have 1 dof per node for the vertical analysis. But the horizontal velocity
    1718                  *runs with 2 dofs per node! Therefore, we multiply by 2 the dof count, to extract the velocity from the inputs: */
    1719                 vx_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+0];
    1720                 vy_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+1];
    1721         }
     1684        if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
     1685            inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
    17221686
    17231687        /*Get gaussian points and weights :*/

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

    r128 r246  
    1212#include "./Node.h"
    1313#include "./Tria.h"
     14#include "./ParameterInputs.h"
    1415
    1516class Penta: public Element{
    … …  
    7172                int   MyRank();
    7273                void  Configure(void* loads,void* nodes,void* materials);
    73                 void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    74                 void  CreateKMatrixDiagnosticHoriz( Mat Kgg, ParameterInputs* inputs, int analysis_type);
    75                 void  CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type);
    76                 void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type);
    77                 void  UpdateFromInputs(ParameterInputs* inputs);
     74                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
     75                void  CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type);
     76                void  CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type);
     77                void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     78                void  UpdateFromInputs(void* inputs);
    7879                void  GetDofList(int* doflist,int* pnumberofdofs);
    79                 Matpar* GetMatPar();
     80                void* GetMatPar();
    8081                int   GetShelf();
    81                 void  GetNodes(Node** nodes);
     82                void  GetNodes(void** nodes);
    8283                int   GetOnBed();
    83                 void  Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
    84                 void  Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type);
    85                 void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type);
    86                 void  GradjB(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type);
    87         double Misfit(double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
     84                void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type);
     85                void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
     86                void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
     87                void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
     88        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type);
    8889               
    8990                void          GetThicknessList(double* thickness_list);
    9091                void          GetBedList(double* bed_list);
    9192                Object* copy();
    92                 Tria*  SpawnTria(int g0, int g1, int g2);
     93                void*  SpawnTria(int g0, int g1, int g2);
    9394
    9495                void  GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4);
    … …  
    102103                void  GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4);
    103104                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3l4);
    104                 void  CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs,int analysis_type);
    105                 void  CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type);
     105                void  CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type);
     106                void  CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type);
    106107                void  GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4);
    107108                void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);

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

    r241 r246  
    266266#define __FUNCT__ "Tria::CreateKMatrix"
    267267
    268 void  Tria::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     268void  Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
    269269
    270270        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    … …  
    284284#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz"
    285285
    286 void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     286void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type){
    287287
    288288
    … …  
    328328       
    329329        /*input parameters for structural analysis (diagnostic): */
    330         double* velocity_param=NULL;
    331         double* oldvelocity_param=NULL;
    332         double  vxvy_list[numgrids][2];
    333         double  oldvxvy_list[numgrids][2];
     330        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     331        double  oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    334332        double  thickness;
     333        int     dofs[2]={0,1};
     334
     335        ParameterInputs* inputs=NULL;
     336
     337        /*recover pointers: */
     338        inputs=(ParameterInputs*)vinputs;
    335339
    336340        /*recover extra inputs from users, at current convergence iteration: */
    337         velocity_param=ParameterInputsRecover(inputs,"velocity");
    338         oldvelocity_param=ParameterInputsRecover(inputs,"old_velocity");
     341        inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     342        inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    339343
    340344        /* Get node coordinates and dof list: */
    … …  
    344348        /* Set Ke_gg to 0: */
    345349        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
    346 
    347         /*Initialize vxvy_list: */
    348         for(i=0;i<numgrids;i++){
    349                 if(velocity_param){
    350                         vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
    351                         vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
    352                 }
    353                 else{
    354                         vxvy_list[i][0]=0;
    355                         vxvy_list[i][1]=0;
    356                 }
    357 
    358                 if(oldvelocity_param){
    359                         oldvxvy_list[i][0]=oldvelocity_param[doflist[i*numberofdofspernode+0]];
    360                         oldvxvy_list[i][1]=oldvelocity_param[doflist[i*numberofdofspernode+1]];
    361                 }
    362                 else{
    363                         oldvxvy_list[i][0]=0;
    364                         oldvxvy_list[i][1]=0;
    365                 }
    366         }
    367350
    368351        #ifdef _DEBUGELEMENTS_
    … …  
    516499#undef __FUNCT__
    517500#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction"
    518 void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     501void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type){
    519502
    520503
    … …  
    554537
    555538        /*input parameters for structural analysis (diagnostic): */
    556         double* velocity_param=NULL;
    557         double  vxvy_list[numgrids][2];
     539        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     540        int     dofs[2]={0,1};
    558541
    559542        /*friction: */
    … …  
    564547        double MOUNTAINKEXPONENT=10;
    565548
    566         /*recover extra inputs from users, at current convergence iteration: */
    567         velocity_param=ParameterInputsRecover(inputs,"velocity");
    568 
     549        ParameterInputs* inputs=NULL;
     550
     551        /*recover pointers: */
     552        inputs=(ParameterInputs*)vinputs;
     553       
    569554        /* Get node coordinates and dof list: */
    570555        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    581566        if (friction_type!=2)throw ErrorException(__FUNCT__," non-viscous friction not supported yet!");
    582567
    583         /*Initialize vxvy_list: */
    584         for(i=0;i<numgrids;i++){
    585                 if(velocity_param){
    586                         vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
    587                         vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
    588                 }
    589                 else{
    590                         vxvy_list[i][0]=0;
    591                         vxvy_list[i][1]=0;
    592                 }
    593         }
     568        /*recover extra inputs from users, at current convergence iteration: */
     569        inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    594570
    595571        /*Build alpha2_list used by drag stiffness matrix*/
    … …  
    691667#undef __FUNCT__
    692668#define __FUNCT__ "Tria::CreatePVector"
    693 void  Tria::CreatePVector(Vec pg,ParameterInputs* inputs,int analysis_type){
     669void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type){
    694670       
    695671        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    … …  
    708684#define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"
    709685
    710 void Tria::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type){
     686void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type){
    711687
    712688        int             i,j;
    … …  
    746722        /*input parameters for structural analysis (diagnostic): */
    747723        double  thickness;
     724
     725        ParameterInputs* inputs=NULL;
     726
     727        /*recover pointers: */
     728        inputs=(ParameterInputs*)vinputs;
    748729
    749730        /* Get node coordinates and dof list: */
    … …  
    871852#undef __FUNCT__
    872853#define __FUNCT__ "Tria::UpdateFromInputs"
    873 void  Tria::UpdateFromInputs(ParameterInputs* inputs){
    874 
    875         int i;
    876         int doflist[MAXDOFSPERNODE*3];
    877         int numberofdofspernode;
    878 
    879         /*element: */
    880         double* h_param=NULL;
    881         double* s_param=NULL;
    882         double* b_param=NULL;
    883         double* k_param=NULL;
    884         double* melting_param=NULL;
    885         double* accumulation_param=NULL;
    886        
    887         /*material: */
    888         double* temperature_average_param=NULL;
     854void  Tria::UpdateFromInputs(void* vinputs){
     855
     856        int     dofs[1]={0};
     857        double  temperature_list[3];
    889858        double  temperature_average;
    890         double* B_param=NULL;
     859        double  B_list[3];
    891860        double  B_average;
    892861
    893 
    894         /*Get dof list: */
    895         GetDofList(&doflist[0],&numberofdofspernode);
     862        ParameterInputs* inputs=NULL;
     863
     864        /*recover pointers: */
     865        inputs=(ParameterInputs*)vinputs;
    896866
    897867        /*Update internal data if inputs holds new values: */
    898         h_param=ParameterInputsRecover(inputs,"thickness");
    899         s_param=ParameterInputsRecover(inputs,"surface");
    900         b_param=ParameterInputsRecover(inputs,"bed");
    901         k_param=ParameterInputsRecover(inputs,"drag");
    902         melting_param=ParameterInputsRecover(inputs,"melting");
    903         accumulation_param=ParameterInputsRecover(inputs,"accumulation");
    904 
    905         for(i=0;i<3;i++){
    906                 if(h_param)h[i]=h_param[doflist[i*numberofdofspernode+0]];
    907                 if(s_param)s[i]=s_param[doflist[i*numberofdofspernode+0]];
    908                 if(b_param)b[i]=b_param[doflist[i*numberofdofspernode+0]];
    909                 if(k_param)k[i]=k_param[doflist[i*numberofdofspernode+0]];
    910                 if(melting_param)melting[i]=melting_param[doflist[i*numberofdofspernode+0]];
    911                 if(accumulation_param)accumulation[i]=accumulation_param[doflist[i*numberofdofspernode+0]];
    912         }
    913 
     868        inputs->Recover("thickness",&h[0],1,dofs,3,(void**)nodes);
     869        inputs->Recover("surface",&s[0],1,dofs,3,(void**)nodes);
     870        inputs->Recover("bed",&b[0],1,dofs,3,(void**)nodes);
     871        inputs->Recover("drag",&k[0],1,dofs,3,(void**)nodes);
     872        inputs->Recover("melting",&melting[0],1,dofs,3,(void**)nodes);
     873        inputs->Recover("accumulation_param",&accumulation[0],1,dofs,3,(void**)nodes);
     874       
    914875        //Update material if necessary
    915         temperature_average_param=ParameterInputsRecover(inputs,"temperature_average");
    916         B_param=ParameterInputsRecover(inputs,"B");
    917         if(temperature_average_param){
    918                 temperature_average-0;
    919                 for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
    920                 temperature_average=temperature_average/3.0;
     876        if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){
     877                temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2])/3.0;
    921878                B_average=Paterson(temperature_average);
    922879                matice->SetB(B_average);
    923880        }
    924 
    925         if(B_param){
    926                 B_average=0;
    927                 for(i=0;i<3;i++) B_average+=B_param[doflist[i*numberofdofspernode+0]];
    928                 B_average=B_average/3.0;
     881       
     882        if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){
     883                B_average=(B_list[0]+B_list[1]+B_list[2])/3.0;
    929884                matice->SetB(B_average);
    930885        }
    … …  
    13111266}
    13121267
    1313 Matpar* Tria::GetMatPar(){
     1268void* Tria::GetMatPar(){
    13141269        return matpar;
    13151270}
    … …  
    13201275
    13211276               
    1322 void  Tria::GetNodes(Node** pnodes){
     1277void  Tria::GetNodes(void** vpnodes){
    13231278        int i;
     1279        Node** pnodes=(Node**)vpnodes;
     1280
    13241281        for(i=0;i<3;i++){
    13251282                pnodes[i]=nodes[i];
    … …  
    13461303
    13471304Object* Tria::copy() {
     1305       
    13481306        return new Tria(*this);
     1307
    13491308}
    13501309
    … …  
    13521311#undef __FUNCT__
    13531312#define __FUNCT__ "Tria::Du"
    1354 void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){
     1313void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
    13551314
    13561315        int i;
    … …  
    14041363        double scaley=0;
    14051364        double scale=0;
    1406         double* fit_param=NULL;
    14071365        double  fit=-1;
     1366
     1367        ParameterInputs* inputs=NULL;
     1368
     1369        /*recover pointers: */
     1370        inputs=(ParameterInputs*)vinputs;
    14081371
    14091372        /* Get node coordinates and dof list: */
    … …  
    14151378
    14161379        /* Recover input data: */
    1417         fit_param=ParameterInputsRecover(inputs,"fit");
    1418         if(fit_param){
    1419                 fit=*fit_param;
    1420         }
    1421         else throw ErrorException(__FUNCT__," missing fit input parameter");
     1380        if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter");
    14221381
    14231382        /*Initialize velocities: */
    … …  
    15701529#undef __FUNCT__
    15711530#define __FUNCT__ "Tria::Gradj"
    1572 void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type,char* control_type){
     1531void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){
    15731532
    15741533        if (strcmp(control_type,"drag")==0){
    … …  
    15831542#undef __FUNCT__
    15841543#define __FUNCT__ "Tria::GradjDrag"
    1585 void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){
     1544void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
    15861545
    15871546
    … …  
    16031562        double adjy_list[numgrids];
    16041563
    1605         double* basal_drag_param=NULL;
    1606         double K_list[numgrids];
    16071564        double drag;
    1608 
    1609         double* bed_param=NULL;
    1610         double bed_list[numgrids];
    1611         double* thickness_param=NULL;
    1612         double thickness_list[numgrids];
     1565        int    dofs[1]={0};
    16131566
    16141567        /* gaussian points: */
    … …  
    16501603        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    16511604
     1605        ParameterInputs* inputs=NULL;
     1606
     1607        /*recover pointers: */
     1608        inputs=(ParameterInputs*)vinputs;
     1609
    16521610        /* Get node coordinates and dof list: */
    16531611        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    16581616
    16591617        /* recover input parameters: */
    1660         basal_drag_param=ParameterInputsRecover(inputs,"drag");
    1661         bed_param=ParameterInputsRecover(inputs,"bed");
    1662         thickness_param=ParameterInputsRecover(inputs,"thickness");
     1618        inputs->Recover("drag",&k[0],1,dofs,numgrids,(void**)nodes);
     1619        inputs->Recover("bed",&b[0],1,dofs,numgrids,(void**)nodes);
     1620        inputs->Recover("thickness",&h[0],1,dofs,numgrids,(void**)nodes);
    16631621
    16641622        /*Initialize parameter lists: */
    … …  
    16701628                adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]];
    16711629                adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]];
    1672                 if(basal_drag_param){
    1673                         K_list[i]=basal_drag_param[doflist[i*numberofdofspernode+0]];
    1674                 }
    1675                 else{
    1676                         K_list[i]=k[i];
    1677                 }
    1678                 if(bed_param){
    1679                         bed_list[i]=bed_param[doflist[i*numberofdofspernode+0]];
    1680                 }
    1681                 else{
    1682                         bed_list[i]=b[i];
    1683                 }
    1684                 if(thickness_param){
    1685                         thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]];
    1686                 }
    1687                 else{
    1688                         thickness_list[i]=h[i];
    1689                 }
    16901630        }
    16911631       
    … …  
    17231663                        friction->rho_ice=matpar->GetRhoIce();
    17241664                        friction->rho_water=matpar->GetRhoWater();
    1725                         friction->K=&K_list[0];
    1726                         friction->bed=&bed_list[0];
    1727                         friction->thickness=&thickness_list[0];
     1665                        friction->K=&k[0];
     1666                        friction->bed=&b[0];
     1667                        friction->thickness=&h[0];
    17281668                        friction->velocities=&vxvy_list[0][0];
    17291669                        friction->p=p;
    … …  
    17461686                /*Recover alpha_complement and k: */
    17471687                GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
    1748                 GetParameterValue(&drag, &K_list[0],gauss_l1l2l3);
     1688                GetParameterValue(&drag, &k[0],gauss_l1l2l3);
    17491689                #ifdef _DEBUG_
    17501690                        printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
    … …  
    17961736#undef __FUNCT__
    17971737#define __FUNCT__ "Tria::GradjB"
    1798 void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){
     1738void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
    17991739
    18001740        int i;
    … …  
    18021742        /* node data: */
    18031743        const int    numgrids=3;
    1804         const int    numdofs=2*numgrids;
    1805         const int    NDOF2=2;
     1744        const int    NDOF1=1;
     1745        const int    NDOF2=1;
     1746        const int    numdofs=NDOF1*numgrids;
    18061747        double       xyz_list[numgrids][3];
    18071748        int          doflist[numdofs];
    … …  
    18141755        double adjx_list[numgrids];
    18151756        double adjy_list[numgrids];
    1816 
    1817         double* thickness_param=NULL;
    1818         double thickness_list[numgrids];
    18191757
    18201758        /* gaussian points: */
    … …  
    18501788        double  lambda,mu;
    18511789        double  thickness;
    1852        
     1790        int     dofs[1]={0};
     1791       
     1792        ParameterInputs* inputs=NULL;
     1793
     1794        /*recover pointers: */
     1795        inputs=(ParameterInputs*)vinputs;
    18531796
    18541797        /* Get node coordinates and dof list: */
    … …  
    18601803
    18611804        /* recover input parameters: */
    1862         thickness_param=ParameterInputsRecover(inputs,"thickness");
     1805        inputs->Recover("thickness",&h[0],1,dofs,numgrids,(void**)nodes);
    18631806
    18641807        /*Initialize parameter lists: */
    … …  
    18701813                adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]];
    18711814                adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]];
    1872                 if(thickness_param){
    1873                         thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]];
    1874                 }
    1875                 else{
    1876                         thickness_list[i]=h[i];
    1877                 }
    18781815        }
    18791816       
    … …  
    19401877#undef __FUNCT__
    19411878#define __FUNCT__ "Tria::Misfit"
    1942 double Tria::Misfit(double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){
     1879double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
    19431880        int i;
    19441881       
    … …  
    19791916        double scalex=1;
    19801917        double scaley=1;
    1981         double* fit_param=NULL;
    19821918        double  fit=-1;
     1919
     1920        ParameterInputs* inputs=NULL;
     1921
     1922        /*recover pointers: */
     1923        inputs=(ParameterInputs*)vinputs;
    19831924
    19841925        /* Get node coordinates and dof list: */
    … …  
    19871928       
    19881929        /* Recover input data: */
    1989         fit_param=ParameterInputsRecover(inputs,"fit");
    1990         if(fit_param){
    1991                 fit=*fit_param;
    1992         }
    1993         else ErrorException(__FUNCT__," missing fit input parameter");
     1930        if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter");
    19941931
    19951932        /*Initialize velocities: */
    … …  
    21022039
    21032040
     2041
     2042#undef __FUNCT__
     2043#define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert"
     2044void  Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* vinputs,int analysis_type){
     2045
     2046        int i,j;
     2047
     2048        /* node data: */
     2049        const int    numgrids=3;
     2050        const int    NDOF1=1;
     2051        const int    numdofs=NDOF1*numgrids;
     2052        double       xyz_list[numgrids][3];
     2053        int          doflist[numdofs];
     2054        int          numberofdofspernode;
     2055       
     2056        /* gaussian points: */
     2057        int     num_gauss,ig;
     2058        double* first_gauss_area_coord  =  NULL;
     2059        double* second_gauss_area_coord =  NULL;
     2060        double* third_gauss_area_coord  =  NULL;
     2061        double* gauss_weights           =  NULL;
     2062        double  gauss_weight;
     2063        double  gauss_l1l2l3[3];
     2064
     2065
     2066        /* matrices: */
     2067        double L[1][3];
     2068        double DL_scalar;
     2069
     2070        double Ke_gg[3][3]; //stiffness matrix
     2071        double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point.
     2072        double Jdet;
     2073
     2074        ParameterInputs* inputs=NULL;
     2075
     2076        /*recover pointers: */
     2077        inputs=(ParameterInputs*)vinputs;
     2078
     2079        /* Get node coordinates and dof list: */
     2080        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     2081        GetDofList(&doflist[0],&numberofdofspernode);
     2082
     2083        /* Set Ke_gg to 0: */
     2084        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
     2085
     2086        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2087        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2088
     2089        /* Start  looping on the number of gaussian points: */
     2090        for (ig=0; ig<num_gauss; ig++){
     2091                /*Pick up the gaussian point: */
     2092                gauss_weight=*(gauss_weights+ig);
     2093                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2094                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2095                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2096
     2097                /* Get Jacobian determinant: */
     2098                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2099       
     2100                //Get L matrix if viscous basal drag present:
     2101                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     2102
     2103                DL_scalar=gauss_weight*Jdet;
     2104               
     2105                /*  Do the triple producte tL*D*L: */
     2106                TripleMultiply( &L[0][0],1,3,1,
     2107                                        &DL_scalar,1,1,0,
     2108                                        &L[0][0],1,3,0,
     2109                                        &Ke_gg_gaussian[0][0],0);
     2110
     2111               
     2112                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     2113                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2114        } //for (ig=0; ig<num_gauss; ig++
     2115
     2116        /*Add Ke_gg to global matrix Kgg: */
     2117        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     2118               
     2119        cleanup_and_return:
     2120        xfree((void**)&first_gauss_area_coord);
     2121        xfree((void**)&second_gauss_area_coord);
     2122        xfree((void**)&third_gauss_area_coord);
     2123        xfree((void**)&gauss_weights);
     2124}
     2125
    21042126#undef __FUNCT__
    21052127#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
    2106 void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     2128void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type){
    21072129
    21082130        int i,j;
    … …  
    21442166        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix
    21452167        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
     2168
     2169        ParameterInputs* inputs=NULL;
     2170
     2171        /*recover pointers: */
     2172        inputs=(ParameterInputs*)vinputs;
    21462173       
    21472174        /* Get node coordinates and dof list: */
    … …  
    22272254#undef __FUNCT__
    22282255#define __FUNCT__ "Tria::CreatePVectorDiagnosticBaseVert"
    2229 void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,ParameterInputs* inputs,int analysis_type){
     2256void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,void* vinputs,int analysis_type){
    22302257
    22312258        int             i,j;
    … …  
    22632290        /*input parameters for structural analysis (diagnostic): */
    22642291        double* velocity_param=NULL;
    2265         double  vx_list[numgrids];
    2266         double  vy_list[numgrids];
     2292        double  vx_list[numgrids]={0,0,0};
     2293        double  vy_list[numgrids]={0,0,0};
    22672294        double  vx,vy;
    22682295        double  meltingvalue;
    22692296        double  slope[2];
    22702297        double  dbdx,dbdy;
     2298        int     dofs1[1]={0};
     2299        int     dofs2[1]={1};
     2300
     2301        ParameterInputs* inputs=NULL;
     2302
     2303        /*recover pointers: */
     2304        inputs=(ParameterInputs*)vinputs;
    22712305
    22722306        /* recover input parameters: */
    2273         velocity_param=ParameterInputsRecover(inputs,"velocity");
    2274         if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity");
     2307        if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity");
     2308            inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
    22752309
    22762310        /* Get node coordinates and dof list: */
    … …  
    22812315        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
    22822316
    2283         /*Initialize vx_list and vy_list: */
    2284         for(i=0;i<numgrids;i++){
    2285                 if(velocity_param){
    2286                         /*doflist is a list of dofs for the node. We have 1 dof per node for the vertical analysis. But the horizontal velocity
    2287                          *runs with 2 dofs per node! Therefore, we multiply by 2 the dof count, to extract the velocity from the inputs: */
    2288                         vx_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+0];
    2289                         vy_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+1];
    2290                 }
    2291                 else{
    2292                         vx_list[i]=0;
    2293                         vy_list[i]=0;
    2294                 }
    2295         }
    2296                
    22972317        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    22982318        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);

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

    r241 r246  
    33 */
    44
    5 #ifndef _TRIA_H
    6 #define _TRIA_H
     5#ifndef _TRIA_H_
     6#define _TRIA_H_
     7
    78
    89#include "./Element.h"
    … …  
    1011#include "./Matice.h"
    1112#include "./Matpar.h"
     13#include "./ParameterInputs.h"
    1214
    1315class Tria: public Element{
    … …  
    6365                int   MyRank();
    6466                void  Configure(void* loads,void* nodes,void* materials);
    65                 void  CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    66                 void  CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type);
    67                 void  UpdateFromInputs(ParameterInputs* inputs);
     67                void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
     68                void  CreatePVector(Vec pg, void* inputs, int analysis_type);
     69                void  UpdateFromInputs(void* inputs);
    6870                void  GetDofList(int* doflist,int* pnumberofdofs);
    6971               
    70                 void  CreateKMatrixDiagnosticHoriz(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    71                 void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,ParameterInputs* inputs,int analysis_type);
    72                 void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,ParameterInputs* inputs,int analysis_type);
     72                void  CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type);
     73                void  CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type);
     74                void  CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* inputs,int analysis_type);
     75                void  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type);
    7376                void  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
    7477                void  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
    … …  
    8386                void  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3);
    8487                void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
    85                 void  Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
    86                 void  Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type);
    87                 void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type);
    88                 void  GradjB(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type);
    89         double Misfit(double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
     88                void  Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type);
     89                void  Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type);
     90                void  GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
     91                void  GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type);
     92        double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type);
    9093
    91                 void  CreatePVectorDiagnosticHoriz(Vec pg,ParameterInputs* inputs,int analysis_type);
    92                 void  CreatePVectorDiagnosticBaseVert(Vec pg,ParameterInputs* inputs,int analysis_type);
    93                 Matpar* GetMatPar();
     94                void  CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type);
     95                void  CreatePVectorDiagnosticBaseVert(Vec pg,void* inputs,int analysis_type);
     96                void* GetMatPar();
    9497                int   GetShelf();
    95                 void  GetNodes(Node** nodes);
     98                void  GetNodes(void** nodes);
    9699                int   GetOnBed();
    97100

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

    r8 r246  
    1414#include "./Matpar.h"
    1515#include "./Node.h"
    16 #include "./Object.h"
    1716#include "./Penta.h"
    1817#include "./Tria.h"
    … …  
    2827#include "./Contour.h"
    2928#include "./ParameterInputs.h"
     29#include "./Input.h"
    3030#include "./Friction.h"
    3131#include "./SolverEnum.h"

  • issm/trunk/src/c/parallel/control.cpp

    r237 r246  
    3434        double search_scalar;
    3535        char* control_type=NULL;
    36         int   gsize;
    3736        double* fit=NULL;
    3837        double* optscal=NULL;
    3938        double* u_g_obs=NULL;
     39        double* u_g_initial=NULL;
     40        int  gsize;
     41        int  numberofnodes;
    4042        double*    maxiter=NULL;
    4143        double  tolx;
    … …  
    5153        OptArgs   optargs;
    5254        OptPars   optpars;
     55        Param*    param=NULL;
    5356
    5457
    … …  
    8790        femmodel.parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint");
    8891        femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint");
    89         gsize=femmodel.nodes->NumberOfDofs();
    90        
     92        femmodel.parameters->FindParam((void*)&gsize,"gsize");
    9193        femmodel.parameters->FindParam((void*)&p_g,"p_g");
    9294        femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs");
    … …  
    9698
    9799        /*Initialize inputs:*/
    98         inputs=NewParameterInputs();
     100       
     101        femmodel.parameters->FindParam((void*)&u_g_initial,"u_g");
     102        femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     103
     104        inputs=new ParameterInputs;
     105        inputs->Add("velocity",u_g_initial,3,numberofnodes);
     106       
     107        /*erase useless parameters: */
     108        param=(Param*)femmodel.parameters->FindParamObject("p_g");
     109        femmodel.parameters->DeleteObject((Object*)param);
     110
     111        param=(Param*)femmodel.parameters->FindParamObject("u_g_obs");
     112        femmodel.parameters->DeleteObject((Object*)param);
     113
     114        param=(Param*)femmodel.parameters->FindParamObject("u_g");
     115        femmodel.parameters->DeleteObject((Object*)param);
     116
    99117
    100118        /*Start looping: */
    … …  
    102120                       
    103121                _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
    104                 ParameterInputsAddFromMat(inputs,p_g,gsize,control_type);
    105                 ParameterInputsAddFromDouble(inputs,fit[n],"fit");
     122                inputs->Add(control_type,p_g,2,numberofnodes);
     123                inputs->Add("fit",fit[n]);
    106124
    107125                _printf_("%s\n","      computing gradJ...");
    … …  
    135153                if (((n+1)%5)==0){
    136154                        _printf_("%s\n","      saving temporary results...");
    137                         ParameterInputsAddFromMat(inputs,p_g,gsize,control_type);
    138                         ParameterInputsAddFromDouble(inputs,fit[n],"fit");
     155                        inputs->Add(control_type,p_g,2,numberofnodes);
     156                        inputs->Add("fit",fit[n]);
    139157                        diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel);
    140158                        OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
    … …  
    150168        /*Write results to disk: */
    151169        _printf_("%s\n","       preparing final velocity solution...");
    152         ParameterInputsAddFromMat(inputs,p_g,gsize,control_type);
    153         ParameterInputsAddFromDouble(inputs,fit[n],"fit");
     170        inputs->Add(control_type,p_g,2,numberofnodes);
     171        inputs->Add("fit",fit[n]);
    154172       
    155173        diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel);

  • issm/trunk/src/c/parallel/diagnostic.cpp

    r210 r246  
    2424        char* lockname=NULL;
    2525        char* analysis_type="diagnostic_horiz";
     26        int   numberofnodes;
    2627
    2728        /*Intermediary: */
    … …  
    3233       
    3334        double* u_g_initial=NULL;
    34         int     gsize;
    3535        Param*  param=NULL;
    3636
    … …  
    6060
    6161        _printf_("initialize inputs:\n");
    62         inputs=NewParameterInputs();
    63         femmodel.parameters->FindParam((void*)&u_g,"u_g");
    64         femmodel.parameters->FindParam((void*)&gsize,"gsize");
    65         ParameterInputsAddFromMat(inputs,u_g_initial,gsize,"velocity");
     62        femmodel.parameters->FindParam((void*)&u_g_initial,"u_g");
     63        femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    6664       
     65        inputs=new ParameterInputs;
     66        inputs->Add("velocity",u_g_initial,3,numberofnodes);
     67
    6768        param=(Param*)femmodel.parameters->FindParamObject("u_g");
    6869        femmodel.parameters->DeleteObject((Object*)param);

  • issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp

    r81 r246  
    3131        int num_unstable_constraints;
    3232        int count;
     33        int numberofnodes;
    3334
    3435        Vec dug;
    … …  
    4748        fem->parameters->FindParam((void*)&connectivity,"connectivity");
    4849        fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
     50        fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    4951        fem->parameters->FindParam((void*)&solver_string,"solverstring");
    5052        fem->parameters->FindParam((void*)&eps_rel,"eps_rel");
    … …  
    6466
    6567                /*Set input parameters: */
    66                 if(old_ug)ParameterInputsAddFromVec(inputs,old_ug,"old_velocity");
    67                 if(ug)ParameterInputsAddFromVec(inputs,ug,"velocity");
     68                if(old_ug)inputs->Add("old_velocity",old_ug,numberofdofspernode,numberofnodes);
     69                if(ug)inputs->Add("velocity",ug,numberofdofspernode,numberofnodes);
    6870
    6971                //save pointer to old velocity
    … …  
    112114                //Deal with penalty loads
    113115                if (debug) _printf_("   penalty constraints\n");
    114                 ParameterInputsAddFromVec(inputs,ug,"velocity");
    115                 ParameterInputsAddFromVec(inputs,old_ug,"old_velocity");
     116                inputs->Add("old_velocity",old_ug,numberofdofspernode,numberofnodes);
     117                inputs->Add("velocity",ug,numberofdofspernode,numberofnodes);
    116118               
    117119                PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type);
    … …  
    157159
    158160                /*Set input parameters: */
    159                 ParameterInputsAddFromVec(inputs,ug,"velocity");
    160                 ParameterInputsAddFromVec(inputs,old_ug,"old_velocity");
     161                inputs->Add("old_velocity",old_ug,numberofdofspernode,numberofnodes);
     162                inputs->Add("velocity",ug,numberofdofspernode,numberofnodes);
    161163       
    162164                kflag=1; pflag=0; //stiffness generation only

  • issm/trunk/src/c/parallel/objectivefunctionC.cpp

    r46 r246  
    4141        Vec     u_g=NULL;
    4242        double* u_g_double=NULL;
     43        int     numberofnodes;
    4344
    4445
    … …  
    5859        femmodel->parameters->FindParam((void*)&fit,"fit");
    5960        femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
     61        femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    6062
    6163        /*First copy p_g so we don't modify it: */
    … …  
    7072
    7173        /*Add new parameter to inputs: */
    72         ParameterInputsAddFromMat(inputs,p_g_copy,gsize,control_type);
     74        inputs->Add(control_type,p_g_copy,2,numberofnodes);
    7375
    7476        //Run diagnostic with updated parameters.
    … …  
    7779
    7880        //Compute misfit for this velocity field.
    79         ParameterInputsAddFromDouble(inputs,fit[n],"fit");
     81        inputs->Add("fit",fit[n]);
    8082        Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials,
    8183                u_g_double,u_g_obs, inputs,analysis_type);

  • issm/trunk/src/m/solutions/cielo/diagnostic.m

    r202 r246  
    1919
    2020        %initialize inputs
    21         if strcmp(md.analysis_type,'diagnostic_stokes'),
    22                 inputs.velocity=m_dh.parameters.u_g;
    23         else
    24                 inputs.velocity=zeros(2*m_dh.parameters.numberofnodes,1);
    25                 inputs.velocity(1:2:end)=m_dh.parameters.u_g(1:3:end);
    26                 inputs.velocity(2:2:end)=m_dh.parameters.u_g(2:3:end);
    27         end
     21        inputs=parameterinputs;
     22        inputs=add(inputs,'velocity',m_dh.parameters.u_g,'doublevec',3,m_dh.parameters.numberofnodes);
    2823
    2924        %Get horizontal solution.

  • issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m

    r202 r246  
    1212        converged=0;
    1313
    14         [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
    15         if velocity_is_present
    16                 soln(count).u_g=velocity_param;
    17         else   
    18                 soln(count).u_g=[];
    19         end
     14        soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node
    2015        soln(count).u_f=[];
    2116
    … …  
    2722               
    2823                %add velocity to inputs
    29                 inputs.velocity=full(soln(count).u_g);
     24                inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
    3025                if count>1,
    31                         inputs.old_velocity=full(soln(count-1).u_g);
     26                        inputs=add(inputs,'old_velocity',soln(count-1).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
    3227                else
    33                         inputs.old_velocity=full(soln(count).u_g);
     28                        inputs=add(inputs,'old_velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
    3429                end
    3530
    3631                %Update inputs in datasets
    3732                [m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
    38 
     33               
    3934                %system matrices
    4035                [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs);
    … …  
    6762
    6863                %Deal with penalty loads
    69                 inputs.velocity=soln(count).u_g;
     64                inputs=add(inputs,'velocity',soln(count).u_g,'vector',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
    7065       
    7166                %penalty constraints
    … …  
    7368
    7469        %   Figure out if convergence is reached.
    75                 if(count>=3 | velocity_is_present),
    76                         dug=soln(count).u_g-soln(count-1).u_g;
    77                         nduinf=norm(dug,inf)*m.parameters.yts;
    78                         ndu=norm(dug,2);
    79                         nu=norm(soln(count-1).u_g,2);
     70                dug=soln(count).u_g-soln(count-1).u_g;
     71                nduinf=norm(dug,inf)*m.parameters.yts;
     72                ndu=norm(dug,2);
     73                nu=norm(soln(count-1).u_g,2);
    8074
    81                         %Relative criterion
    82                         if (ndu/nu<=m.parameters.eps_rel),
    83                                 if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: norm(du)/norm(u)=',ndu/nu,' < ',m.parameters.eps_rel)); end
    84                                 converged=1;
     75                %Relative criterion
     76                if (ndu/nu<=m.parameters.eps_rel),
     77                        if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: norm(du)/norm(u)=',ndu/nu,' < ',m.parameters.eps_rel)); end
     78                        converged=1;
     79                else
     80                        if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: norm(du)/norm(u)=',ndu/nu,' > ',m.parameters.eps_rel)); end
     81                        converged=0;
     82                end
     83
     84                %Absolute criterion
     85                if ~isnan(m.parameters.eps_abs),
     86                        if (nduinf<=m.parameters.eps_abs),
     87                                if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: max(du)=',nduinf,' < ',m.parameters.eps_abs)); end
    8588                        else
    86                                 if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: norm(du)/norm(u)=',ndu/nu,' > ',m.parameters.eps_rel)); end
    87                                 converged=0;
     89                                if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: max(du)=',nduinf,' > ',m.parameters.eps_abs)); end
     90                                        converged=0;
    8891                        end
     92                end
    8993
    90                         %Absolute criterion
    91                         if ~isnan(m.parameters.eps_abs),
    92                                 if (nduinf<=m.parameters.eps_abs),
    93                                         if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: max(du)=',nduinf,' < ',m.parameters.eps_abs)); end
    94                                 else
    95                                         if m.parameters.debug, disp(sprintf('%s %g %s %g','      convergence criterion: max(du)=',nduinf,' > ',m.parameters.eps_abs)); end
    96                                                 converged=0;
    97                                 end
    98                         end
    99 
    100                         %rift convergence criterion
    101                         if ~constraints_converged,
    102                                 converged=0;
    103                         end
     94                %rift convergence criterion
     95                if ~constraints_converged,
     96                        converged=0;
    10497                end
    10598
    … …  
    118111        if nout==2,
    119112               
    120                 inputs.velocity=full(soln(count).u_g);
     113                inputs=add(inputs,'velocity',soln(count).u_g,'vector',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
    121114                m.parameters.kflag=1; m.parameters.pflag=0;
    122115                [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs);

  • issm/trunk/src/mex/Du/Du.cpp

    r1 r246  
    4040
    4141        /*Fetch inputs: */
    42         ParameterInputsInit(&inputs,INPUTS);
     42        inputs=new ParameterInputs;
     43        inputs->Init(INPUTS);
    4344
    4445        /*!Call core code: */
    … …  
    5657        xfree((void**)&u_g_obs);
    5758        VecFree(&du_g);
    58         DeleteParameterInputs(&inputs);
     59        delete inputs;
    5960
    6061        /*end module: */

  • issm/trunk/src/mex/Gradj/Gradj.cpp

    r1 r246  
    4242
    4343        /*Fetch inputs: */
    44         ParameterInputsInit(&inputs,INPUTS);
     44        inputs=new ParameterInputs;
     45        inputs->Init(INPUTS);
    4546
    4647        /*!Call core code: */
    … …  
    5859        xfree((void**)&lambda_g);
    5960        xfree((void**)&control_type);
    60         DeleteParameterInputs(&inputs);
     61        delete inputs;
    6162        VecFree(&grad_g);
    6263

  • issm/trunk/src/mex/Misfit/Misfit.cpp

    r1 r246  
    3939
    4040        /*Fetch inputs: */
    41         ParameterInputsInit(&inputs,INPUTS);
     41        inputs=new ParameterInputs;
     42        inputs->Init(INPUTS);
    4243
    4344        /*!Call core code: */
    … …  
    5455        xfree((void**)&u_g);
    5556        xfree((void**)&u_g_obs);
    56         DeleteParameterInputs(&inputs);
     57        delete inputs
    5758
    5859        /*end module: */

  • issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp

    r1 r246  
    3838
    3939        /*Fetch inputs: */
    40         ParameterInputsInit(&inputs,INPUTS);
     40        inputs=new ParameterInputs;
     41        inputs->Init(INPUTS);
    4142
    4243        /*!Generate internal degree of freedom numbers: */
    … …  
    5455        delete loads;
    5556        delete materials;
    56         DeleteParameterInputs(&inputs);
     57        delete inputs;
    5758
    5859        /*end module: */

  • issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp

    r1 r246  
    4040
    4141        /*Fetch inputs: */
    42         ParameterInputsInit(&inputs,INPUTS);
     42        inputs=new ParameterInputs;
     43        inputs->Init(INPUTS);
    4344
    4445        /*!Generate stiffnesses from penalties: */
    … …  
    5455        delete loads;
    5556        delete materials;
    56         DeleteParameterInputs(&inputs);
     57        delete inputs;
    5758        MatFree(&Kgg);
    5859        VecFree(&pg);

  • issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp

    r1 r246  
    4444
    4545        /*Fetch inputs: */
    46         ParameterInputsInit(&inputs,INPUTS);
     46        inputs=new ParameterInputs;
     47        inputs->Init(INPUTS);
    4748
    4849        /*!Generate internal degree of freedom numbers: */
    … …  
    5960        delete loads;
    6061        delete materials;
    61         DeleteParameterInputs(&inputs);
     62        delete inputs;
    6263        MatFree(&Kgg);
    6364        VecFree(&pg);

  • issm/trunk/src/mex/UpdateFromInputs/UpdateFromInputs.cpp

    r1 r246  
    3030
    3131        /*Fetch inputs: */
    32         ParameterInputsInit(&inputs,INPUTS);
     32        inputs=new ParameterInputs;
     33        inputs->Init(INPUTS);
    3334
    3435        /*!Generate internal degree of freedom numbers: */
    … …  
    4647        delete loads;
    4748        delete materials;
    48         DeleteParameterInputs(&inputs);
     49        delete inputs;
    4950
    5051        /*end module: */
Note: See TracChangeset for help on using the changeset viewer.