Changeset 202


Ignore:
Timestamp:
05/04/09 10:34:57 (16 years ago)
Author:
Eric.Larour
Message:

Elements types now in setelementstype.m routine. New partitioning (independent of number of dofs

Location:
issm/trunk/src
Files:
1 deleted
16 edited

Legend:

Unmodified
Added
Removed

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

    r100 r202  
    146146}
    147147
    148 int  DataSet::DeleteObject(int id){
    149 
    150         return 0;
     148int  DataSet::DeleteObject(Object* object){
     149       
     150        vector<Object*>::iterator iterator;
     151
     152        iterator = find(objects.begin(), objects.end(),object);
     153
     154        objects.erase(iterator);
     155
    151156}
    152157
    … …  
    293298}
    294299
     300Object*   DataSet::FindParamObject(char* name){
     301
     302        /*Go through a dataset, and find a Param* object
     303         *which parameter name is "name" : */
     304       
     305        vector<Object*>::iterator object;
     306        Param* param=NULL;
     307
     308        for ( object=objects.begin() ; object < objects.end(); object++ ){
     309
     310                /*Find param type objects: */
     311                if((*object)->Enum()==ParamEnum()){
     312
     313                        /*Ok, this object is a parameter, recover it and ask which name it has: */
     314                        param=(Param*)(*object);
     315
     316                        if (strcmp(param->GetParameterName(),name)==0){
     317                                /*Ok, this is the one! Return the object: */
     318                                return (*object);
     319                        }
     320                }
     321        }
     322        return NULL;
     323}
     324
     325
    295326void   DataSet::NodeRank(int* ranks){
    296327
    … …  
    324355
    325356        int  dofcount=0;
     357        int  dofcount1=0;
    326358        int* alldofcount=NULL;
     359        int* alldofcount1=NULL;
    327360        int* borderdofs=NULL;
     361        int* borderdofs1=NULL;
    328362        int* allborderdofs=NULL;
     363        int* allborderdofs1=NULL;
    329364        int  i;
    330365        vector<Object*>::iterator object;
    … …  
    343378                       
    344379                        /*Ok, this object is a node, ask it to distribute dofs, and update the dofcount: */
    345                         node->DistributeDofs(&dofcount);
     380                        node->DistributeDofs(&dofcount,&dofcount1);
    346381                       
    347382                }
    … …  
    356391        MPI_Bcast(alldofcount,num_procs,MPI_INT,0,MPI_COMM_WORLD);
    357392
     393        alldofcount1=(int*)xmalloc(num_procs*sizeof(int));
     394        MPI_Gather(&dofcount1,1,MPI_INT,alldofcount1,1,MPI_INT,0,MPI_COMM_WORLD);
     395        MPI_Bcast(alldofcount1,num_procs,MPI_INT,0,MPI_COMM_WORLD);
     396
    358397        /*Ok, now every cpu should start its own dof count at the end of the dofcount
    359398         * from cpu-1. : */
    360399        dofcount=0;
     400        dofcount1=0;
    361401        if(my_rank==0){
    362402                dofcount=0;
     403                dofcount1=0;
    363404        }
    364405        else{
    365406                for(i=0;i<my_rank;i++){
    366407                        dofcount+=alldofcount[i];
     408                        dofcount1+=alldofcount1[i];
    367409                }
    368410        }
    … …  
    378420                       
    379421                        /*Ok, this object is a node, ask it to update his dofs: */
    380                         node->UpdateDofs(dofcount);
     422                        node->UpdateDofs(dofcount,dofcount1);
    381423                       
    382424                }
    … …  
    387429        borderdofs=(int*)xcalloc(numberofnodes*numdofspernode,sizeof(int));
    388430        allborderdofs=(int*)xcalloc(numberofnodes*numdofspernode,sizeof(int));
     431        borderdofs1=(int*)xcalloc(numberofnodes*3,sizeof(int));
     432        allborderdofs1=(int*)xcalloc(numberofnodes*3,sizeof(int));
     433
    389434        for ( object=objects.begin() ; object < objects.end(); object++ ){
    390435
    … …  
    395440                       
    396441                        /*Ok, let this object show its border dofs, if is is a border dof: */
    397                         node->ShowBorderDofs(borderdofs);
     442                        node->ShowBorderDofs(borderdofs,borderdofs1);
    398443                       
    399444                }
    400445        }
    401446        MPI_Allreduce ( (void*)borderdofs,(void*)allborderdofs,numberofnodes*numdofspernode,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
     447        MPI_Allreduce ( (void*)borderdofs1,(void*)allborderdofs1,numberofnodes*3,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
    402448
    403449        /*Ok, now every cpu knows everyone else's border node dofs, update the border dofs accordingly: */
    … …  
    410456                       
    411457                        /*Ok, this object is a node, ask it to update his dofs: */
    412                         node->UpdateBorderDofs(allborderdofs);
     458                        node->UpdateBorderDofs(allborderdofs,allborderdofs1);
    413459                       
    414460                }
    … …  
    420466        xfree((void**)&borderdofs);
    421467        xfree((void**)&allborderdofs);
     468        xfree((void**)&alldofcount1);
     469        xfree((void**)&borderdofs1);
     470        xfree((void**)&allborderdofs1);
     471
    422472        return;
    423473
    … …  
    433483
    434484        /*Create partition vector: */
    435         partition=NewVec(numberofnodes*numdofspernode);
    436 
    437         /*Go through all nodes, and ask each node to plug its doflist in
     485        partition=NewVec(numberofnodes);
     486
     487        /*Go through all nodes, and ask each node to plug its 1D doflist in
    438488         * partition. The location where each node plugs its doflist into
    439          * partition is determined by its (id-1)*numdofspernode (ie, serial
    440          * organisation of the dofs).
     489         * partition is determined by its (id-1)*3 (ie, serial * organisation of the dofs).
    441490         */
    442491       

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

    r100 r202  
    4545                int   Size();
    4646                int   FindParam(void* pvalue, char* name);
     47                Object* FindParamObject(char* name);
    4748                void  NodeRank(int* ranks);
    4849                void  DistributeDofs(int numberofnodes,int numdofspernode);
    … …  
    7677                void  Misfit(double* pJ, double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type);
    7778                void  VelocityExtrude(Vec ug,double* ug_serial);
     79                int   DeleteObject(Object* object);
    7880
    7981

  • issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp

    r117 r202  
    2525        double* optscal=NULL;
    2626        double* maxiter=NULL;
    27         double* control_parameter=NULL;
     27        Vec     control_parameter=NULL;
     28        Vec     vx_obs=NULL;
     29        Vec     vy_obs=NULL;
    2830
    2931        /*Get counter : */
    … …  
    98100
    99101        /*Get vx_obs, vy_obs, and the parameter value: */
    100         ModelFetchData((void**)&model->vx_obs,NULL,NULL,model_handle,"vx_obs","Matrix","Mat");
    101         ModelFetchData((void**)&model->vy_obs,NULL,NULL,model_handle,"vy_obs","Matrix","Mat");
    102         ModelFetchData((void**)&control_parameter,NULL,NULL,model_handle,model->control_type,"Matrix","Mat");
     102        ModelFetchData((void**)&vx_obs,NULL,NULL,model_handle,"vx_obs","Vector",NULL);
     103        ModelFetchData((void**)&vy_obs,NULL,NULL,model_handle,"vy_obs","Vector",NULL);
     104        ModelFetchData((void**)&control_parameter,NULL,NULL,model_handle,model->control_type,"Vector",NULL);
    103105
    104         for(i=0;i<model->numberofnodes;i++)model->vx_obs[i]=model->vx_obs[i]/model->yts;
    105         for(i=0;i<model->numberofnodes;i++)model->vy_obs[i]=model->vy_obs[i]/model->yts;
     106        VecScale(vx_obs,1.0/model->yts);
     107        VecScale(vy_obs,1.0/model->yts);
    106108
    107109        count++;
    108         param= new Param(count,"vx_obs",DOUBLEVEC);
    109         param->SetDoubleVec(model->vx_obs,model->numberofnodes);
     110        param= new Param(count,"vx_obs",PETSCVEC);
     111        param->SetVec(vx_obs);
    110112        parameters->AddObject(param);
    111113
    112114        count++;
    113         param= new Param(count,"vy_obs",DOUBLEVEC);
    114         param->SetDoubleVec(model->vy_obs,model->numberofnodes);
     115        param= new Param(count,"vy_obs",PETSCVEC);
     116        param->SetVec(vy_obs);
    115117        parameters->AddObject(param);
    116118
    117119        count++;
    118         param= new Param(count,"control_parameter",DOUBLEVEC);
    119         param->SetDoubleVec(control_parameter,model->numberofnodes);
     120        param= new Param(count,"control_parameter",PETSCVEC);
     121        param->SetVec(control_parameter);
    120122        parameters->AddObject(param);
    121123       
    122         xfree((void**)&model->vx_obs);
    123         xfree((void**)&model->vy_obs);
    124         xfree((void**)&control_parameter);
     124        VecFree(&vx_obs);
     125        VecFree(&vy_obs);
     126        VecFree(&control_parameter);
    125127
    126128

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

    r117 r202  
    2727        double* maxiter=NULL;
    2828        double* control_parameter=NULL;
    29        
    30        
     29        Vec     vx=NULL;
     30        Vec     vy=NULL;
     31        Vec     vz=NULL;
    3132       
    3233        /*Initialize dataset: */
    … …  
    168169        param->SetInteger(numberofdofspernode);
    169170        parameters->AddObject(param);
    170 
    171171
    172172        /*Deal with control parameters: */
    … …  
    175175        }
    176176
     177        /*Diagnostic: */
     178        /*Get vx and vy: */
     179        ModelFetchData((void**)&vx,NULL,NULL,model_handle,"vx","Vector",NULL);
     180        ModelFetchData((void**)&vy,NULL,NULL,model_handle,"vy","Vector",NULL);
     181        ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Vector",NULL);
     182
     183        if(vx)VecScale(vx,1.0/model->yts);
     184        if(vy)VecScale(vy,1.0/model->yts);
     185        if(vz)VecScale(vz,1.0/model->yts);
     186
     187        count++;
     188        param= new Param(count,"vx",PETSCVEC);
     189        param->SetVec(vx);
     190        parameters->AddObject(param);
     191
     192        count++;
     193        param= new Param(count,"vy",PETSCVEC);
     194        param->SetVec(vy);
     195        parameters->AddObject(param);
     196
     197        count++;
     198        param= new Param(count,"vz",PETSCVEC);
     199        param->SetVec(vz);
     200        parameters->AddObject(param);
     201
     202        VecFree(&vx);
     203        VecFree(&vy);
     204        VecFree(&vz);
     205
    177206        /*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these
    178207         * datasets, it will not be redone: */

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

    r59 r202  
    1919        double* partition=NULL;
    2020        Param*  param=NULL;
     21
     22        /*diagnostic: */
     23        Vec     vec_vx=NULL;
     24        Vec     vec_vy=NULL;
     25        Vec     vec_vz=NULL;
     26
     27        double* vx=NULL;
     28        double* vy=NULL;
     29        double* vz=NULL;
     30       
     31        double* u_g=NULL;
     32       
     33        /*control: */
     34        Vec     vec_vx_obs=NULL;
     35        Vec     vec_vy_obs=NULL;
     36        Vec     vec_control_parameter=NULL;
     37
    2138        double* vx_obs=NULL;
    2239        double* vy_obs=NULL;
    2340        double* control_parameter=NULL;
     41
     42        double* u_g_obs=NULL;
     43        double* p_g=NULL;
     44       
     45
     46        /*diverse: */
    2447        int     numberofnodes;
    2548        int     analysis_type;
    2649        int     count;
    2750
    28         /*new parameters: */
    29         double* u_g_obs=NULL;
    30         double* p_g=NULL;
    31 
    3251        parameters->FindParam((void*)&analysis_type,"analysis_type");
    3352        count=parameters->Size();
     53               
     54        parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     55
     56        /*First serialize partition vector: */
     57        VecToMPISerial(&partition,part);
     58       
     59       
     60        if (   (analysis_type==ControlAnalysisEnum()) ||
     61               (analysis_type==DiagnosticHorizAnalysisEnum()) ||
     62               (analysis_type==DiagnosticVertAnalysisEnum()) ||
     63               (analysis_type==DiagnosticStokesAnalysisEnum())
     64           ){
     65
     66                parameters->FindParam((void*)&vec_vx,"vx");
     67                parameters->FindParam((void*)&vec_vy,"vy");
     68                parameters->FindParam((void*)&vec_vz,"vz");
     69
     70                if(vec_vx)VecToMPISerial(&vx,vec_vx);
     71                if(vec_vy)VecToMPISerial(&vy,vec_vy);
     72                if(vec_vz)VecToMPISerial(&vz,vec_vz);
     73                       
     74                u_g=(double*)xcalloc(numberofnodes*3,sizeof(double));
     75
     76                if(vx){
     77                        for(i=0;i<numberofnodes;i++){
     78                                u_g[(int)(3*partition[i]+0)]=vx[i]; 
     79                        }
     80                }
     81
     82                if(vy){
     83                        for(i=0;i<numberofnodes;i++){
     84                                u_g[(int)(3*partition[i]+1)]=vy[i]; 
     85                        }
     86                }
     87
     88                if(vz){
     89                        for(i=0;i<numberofnodes;i++){
     90                                u_g[(int)(3*partition[i]+2)]=vz[i]; 
     91                        }
     92                }
     93
     94
     95                /*Now, create new parameters: */
     96                count++;
     97                param= new Param(count,"u_g",DOUBLEVEC);
     98                param->SetDoubleVec(u_g,3*numberofnodes);
     99                parameters->AddObject(param);
     100
     101                /*erase old parameters: */
     102                param=(Param*)parameters->FindParamObject("vx");
     103                parameters->DeleteObject((Object*)param);
     104
     105                param=(Param*)parameters->FindParamObject("vy");
     106                parameters->DeleteObject((Object*)param);
     107
     108                param=(Param*)parameters->FindParamObject("vz");
     109                parameters->DeleteObject((Object*)param);
     110
     111        }
     112
    34113
    35114        if(analysis_type==ControlAnalysisEnum()){
    36115
    37                 /*First serialize partition vector: */
    38                 VecToMPISerial(&partition,part);
    39                
    40                 parameters->FindParam((void*)&vx_obs,"vx_obs");
    41                 parameters->FindParam((void*)&vy_obs,"vy_obs");
    42                 parameters->FindParam((void*)&control_parameter,"control_parameter");
    43                 parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     116                parameters->FindParam((void*)&vec_vx_obs,"vx_obs");
     117                parameters->FindParam((void*)&vec_vy_obs,"vy_obs");
     118                parameters->FindParam((void*)&vec_control_parameter,"control_parameter");
    44119
    45120                /*Now, from vx_obs and vy_obs, build u_g_obs, correctly partitioned: */
     121                VecToMPISerial(&vx_obs,vec_vx_obs);
     122                VecToMPISerial(&vy_obs,vec_vy_obs);
     123                VecToMPISerial(&control_parameter,vec_control_parameter);
     124
    46125                u_g_obs=(double*)xcalloc(numberofnodes*2,sizeof(double));
    47126                p_g=(double*)xcalloc(numberofnodes*2,sizeof(double));
    48127
    49128                for(i=0;i<numberofnodes;i++){
    50                         u_g_obs[(int)(partition[2*i+0])]=vx_obs[i]; 
    51                         u_g_obs[(int)(partition[2*i+1])]=vy_obs[i]; 
    52                         p_g[(int)(partition[2*i+0])]=control_parameter[i];
     129                        u_g_obs[(int)(2*partition[i]+0)]=vx_obs[i]; 
     130                        u_g_obs[(int)(2*partition[i]+1)]=vy_obs[i]; 
     131                        p_g[(int)(2*partition[i]+0)]=control_parameter[i];
    53132                }
    54133
    … …  
    61140                count++;
    62141                param= new Param(count,"p_g",DOUBLEVEC);
    63                 param->SetDoubleVec(p_g,2*numberofnodes);
     142                param->SetDoubleVec(p_g,numberofnodes);
    64143                parameters->AddObject(param);
     144
     145                /*erase old parameter: */
     146                param=(Param*)parameters->FindParamObject("vx_obs");
     147                parameters->DeleteObject((Object*)param);
     148
     149                param=(Param*)parameters->FindParamObject("vy_obs");
     150                parameters->DeleteObject((Object*)param);
     151
     152                param=(Param*)parameters->FindParamObject("control_parameter");
     153                parameters->DeleteObject((Object*)param);
     154
    65155        }
    66156
    67157        xfree((void**)&partition);
     158       
     159        xfree((void**)&vx);
     160        xfree((void**)&vy);
     161        xfree((void**)&vz);
     162        VecFree(&vec_vx);
     163        VecFree(&vec_vy);
     164        VecFree(&vec_vz);
     165        xfree((void**)&u_g);
     166
    68167        xfree((void**)&vx_obs);
    69168        xfree((void**)&vy_obs);
     169        VecFree(&vec_vx_obs);
     170        VecFree(&vec_vy_obs);
     171        VecFree(&vec_control_parameter);
    70172        xfree((void**)&control_parameter);
    71173        xfree((void**)&u_g_obs);
    72174        xfree((void**)&p_g);
     175
    73176}

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

    r137 r202  
    7373                printf("%i|",doflist[i]);
    7474        }
     75        printf("   doflist1:|");
     76        printf("%i|",doflist1[1]);
    7577        printf("\n");
    7678
    … …  
    109111        memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
    110112        memcpy(marshalled_dataset,&doflist,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
     113        memcpy(marshalled_dataset,&doflist1,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
    111114        memcpy(marshalled_dataset,&mset,sizeof(mset));marshalled_dataset+=sizeof(mset);
    112115        memcpy(marshalled_dataset,&nset,sizeof(nset));marshalled_dataset+=sizeof(nset);
    … …  
    131134                sizeof(onsurface)+
    132135                sizeof(doflist)+
     136                sizeof(doflist1)+
    133137                sizeof(mset)+
    134138                sizeof(nset)+
    … …  
    164168        memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
    165169        memcpy(&doflist,marshalled_dataset,sizeof(doflist));marshalled_dataset+=sizeof(doflist);
     170        memcpy(&doflist1,marshalled_dataset,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1);
    166171        memcpy(&mset,marshalled_dataset,sizeof(mset));marshalled_dataset+=sizeof(mset);
    167172        memcpy(&nset,marshalled_dataset,sizeof(nset));marshalled_dataset+=sizeof(nset);
    … …  
    193198}
    194199
    195 void  Node::DistributeDofs(int* pdofcount){
     200void  Node::DistributeDofs(int* pdofcount,int* pdofcount1){
    196201
    197202        int i;
    198203        extern int my_rank;
    199204        int dofcount;
     205        int dofcount1;
    200206
    201207        dofcount=*pdofcount;
     208        dofcount1=*pdofcount1;
    202209       
    203210        if(clone){
    … …  
    212219        dofcount+=numberofdofs;
    213220
     221        doflist1[0]=dofcount1;
     222        dofcount1+=1;
     223
    214224        /*Assign output pointers: */
    215225        *pdofcount=dofcount;
     226        *pdofcount1=dofcount1;
    216227
    217228        return;
    218229}
    219230
    220 void  Node::UpdateDofs(int dofcount){
     231void  Node::UpdateDofs(int dofcount,int dofcount1){
    221232       
    222233        int i;
    … …  
    232243                doflist[i]+=dofcount;
    233244        }
     245        doflist1[0]+=dofcount1;
    234246
    235247        return;
    236248}
    237249
    238 void  Node::ShowBorderDofs(int* borderdofs){
     250void  Node::ShowBorderDofs(int* borderdofs,int* borderdofs1){
    239251
    240252        int j;
    … …  
    252264                *(borderdofs+numberofdofs*(id-1)+j)=doflist[j];
    253265        }
    254 
    255         return;
    256 }
    257 
    258 void  Node::UpdateBorderDofs(int* allborderdofs){
     266        *(borderdofs1+(id-1)+0)=doflist1[0];
     267
     268        return;
     269}
     270
     271void  Node::UpdateBorderDofs(int* allborderdofs,int* allborderdofs1){
    259272
    260273        int j;
    … …  
    272285        for(j=0;j<numberofdofs;j++){
    273286                doflist[j]=*(allborderdofs+numberofdofs*(id-1)+j);
    274         }       
     287        }
     288        doflist1[0]=*(allborderdofs1+(id-1)+0);
    275289        return;
    276290}
    277291void  Node::CreatePartition(Vec partition){
    278292
    279         int      i;
    280         int*     idxm=NULL;
    281         double*  values=NULL;
    282 
    283         idxm=(int*)xmalloc(numberofdofs*sizeof(int));
    284         values=(double*)xmalloc(numberofdofs*sizeof(double));
    285 
    286         for(i=0;i<numberofdofs;i++)idxm[i]=(id-1)*numberofdofs+i;
    287         for(i=0;i<numberofdofs;i++)values[i]=(double)doflist[i];
    288 
    289         VecSetValues(partition,numberofdofs,idxm,values,INSERT_VALUES);
    290 
    291         xfree((void**)&idxm);
    292         xfree((void**)&values);
    293 
    294         return;
    295 }
    296 void  Node::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
    297                
     293        int      idxm;
     294        double   value;
     295
     296        idxm=(id-1);
     297        value=(double)doflist1[0];
     298
     299        VecSetValues(partition,1,&idxm,&value,INSERT_VALUES);
     300
     301        return;
     302}
    298303               
    299304void  Node::SetClone(int* minranks){

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

    r95 r202  
    3636
    3737                /*data that is post processed : */
    38                 int doflist[MAXDOFSPERNODE];
     38                int doflist[MAXDOFSPERNODE]; //dof list on which we solve
     39                int doflist1[1]; //1d dof list to recover input parameters
    3940
    4041        public:
    … …  
    5253                int   GetId(void);
    5354                int   MyRank(void);
    54                 void  DistributeDofs(int* pdofcount);
    55                 void  UpdateDofs(int dofcount);
    56                 void  ShowBorderDofs(int* borderdofs);
    57                 void  UpdateBorderDofs(int* allborderdofs);
     55                void  DistributeDofs(int* pdofcount,int* pdofcount1);
     56                void  UpdateDofs(int dofcount,int dofcount1);
     57                void  ShowBorderDofs(int* borderdofs,int* borderdofs1);
     58                void  UpdateBorderDofs(int* allborderdofs,int* allborderdofs1);
    5859                void  CreatePartition(Vec partition);
    59                 void  DistributeNumDofs(int* numdofspernode,int analysis_type);
    6060                void  SetClone(int* minranks);
    6161                int   GetNumberOfDofs();

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

    r1 r202  
    3737
    3838Param::~Param(){
    39         return;
     39
     40        switch(type){
     41       
     42                case STRING:
     43                        break;
     44       
     45                case INTEGER:
     46                        break;
     47       
     48                case DOUBLE:
     49                        break;
     50               
     51                case DOUBLEVEC:
     52                        xfree((void**)&doublevec);
     53                        break;
     54       
     55                case DOUBLEMAT:
     56                        xfree((void**)&doublemat);
     57                        break;
     58
     59                case PETSCVEC:
     60                        VecFree(&vec);
     61                        break;
     62
     63                case  PETSCMAT:
     64                        MatFree(&mat);
     65                        break;
     66
     67                default:
     68                        throw ErrorException(__FUNCT__,exprintf("%s%i","unknow parameter type ",type));
     69        }
    4070}
    4171
    … …  
    139169
    140170                case PETSCVEC:
    141                         VecGetSize(vec,&M);
    142                         VecToMPISerial(&serial_vec,vec);
    143                         memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
    144                         memcpy(marshalled_dataset,serial_vec,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
    145                         xfree((void**)&serial_vec);
     171                        if(vec){
     172                                VecGetSize(vec,&M);
     173                                VecToMPISerial(&serial_vec,vec);
     174                                memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
     175                                memcpy(marshalled_dataset,serial_vec,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
     176                                xfree((void**)&serial_vec);
     177                        }
     178                        else{
     179                                M=0;
     180                                memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
     181                        }
    146182                        break;
    147183
    148184                case PETSCMAT:
    149                         MatGetSize(mat,&M,&N);
    150                         MatToSerial(&serial_mat,mat);
    151                         memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
    152                         memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
    153                         memcpy(marshalled_dataset,serial_mat,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
    154                         xfree((void**)&serial_mat);
     185                        if(mat){
     186                                MatGetSize(mat,&M,&N);
     187                                MatToSerial(&serial_mat,mat);
     188                                memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
     189                                memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
     190                                memcpy(marshalled_dataset,serial_mat,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
     191                                xfree((void**)&serial_mat);
     192                        }
     193                        else{
     194                                M=0;
     195                                N=0;
     196                                memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
     197                                memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
     198                        }
    155199                        break;
    156200               
    … …  
    263307                case PETSCVEC:
    264308                        memcpy(&M,marshalled_dataset,sizeof(M));marshalled_dataset+=sizeof(M);
    265                         serial_vec=(double*)xmalloc(M*sizeof(double));
    266                         memcpy(serial_vec,marshalled_dataset,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
    267 
    268                         vec=NewVec(M);
    269                         idxm=(int*)xmalloc(M*sizeof(int));
    270                         for(i=0;i<M;i++)idxm[i]=i;
    271                         VecSetValues(vec,M,idxm,serial_vec,INSERT_VALUES);
    272 
    273                         VecAssemblyBegin(vec);
    274                         VecAssemblyEnd(vec);
    275                        
    276                         xfree((void**)&serial_vec);
    277                         xfree((void**)&idxm);
     309                        if(M){
     310                                serial_vec=(double*)xmalloc(M*sizeof(double));
     311                                memcpy(serial_vec,marshalled_dataset,M*sizeof(double));marshalled_dataset+=(M*sizeof(double));
     312
     313                                vec=NewVec(M);
     314                                idxm=(int*)xmalloc(M*sizeof(int));
     315                                for(i=0;i<M;i++)idxm[i]=i;
     316                                VecSetValues(vec,M,idxm,serial_vec,INSERT_VALUES);
     317
     318                                VecAssemblyBegin(vec);
     319                                VecAssemblyEnd(vec);
     320                               
     321                                xfree((void**)&serial_vec);
     322                                xfree((void**)&idxm);
     323                        }
     324                        else{
     325                                vec=NULL;
     326                        }
    278327                        break;
    279328
    … …  
    281330                        memcpy(&M,marshalled_dataset,sizeof(M));marshalled_dataset+=sizeof(M);
    282331                        memcpy(&N,marshalled_dataset,sizeof(N));marshalled_dataset+=sizeof(N);
    283                         serial_mat=(double*)xmalloc(M*N*sizeof(double));
    284                         memcpy(serial_mat,marshalled_dataset,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
    285                         mat=NewMat(M,N,&sparsity,NULL,NULL);
    286                         idxm=(int*)xmalloc(M*sizeof(int));
    287                         idxn=(int*)xmalloc(N*sizeof(int));
    288                         for(i=0;i<M;i++)idxm[i]=i;
    289                         for(i=0;i<N;i++)idxn[i]=i;
    290                         MatSetValues(mat,M,idxm,N,idxn,serial_mat,INSERT_VALUES);
    291                         MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
    292                         MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
    293 
    294                         xfree((void**)&serial_mat);
    295                         xfree((void**)&idxm);
    296                         xfree((void**)&idxn);
     332                        if(M!=0 && N!=0){
     333                                serial_mat=(double*)xmalloc(M*N*sizeof(double));
     334                                memcpy(serial_mat,marshalled_dataset,M*N*sizeof(double));marshalled_dataset+=(M*N*sizeof(double));
     335                                mat=NewMat(M,N,&sparsity,NULL,NULL);
     336                                idxm=(int*)xmalloc(M*sizeof(int));
     337                                idxn=(int*)xmalloc(N*sizeof(int));
     338                                for(i=0;i<M;i++)idxm[i]=i;
     339                                for(i=0;i<N;i++)idxn[i]=i;
     340                                MatSetValues(mat,M,idxm,N,idxn,serial_mat,INSERT_VALUES);
     341                                MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
     342                                MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
     343
     344                                xfree((void**)&serial_mat);
     345                                xfree((void**)&idxm);
     346                                xfree((void**)&idxn);
     347                        }
     348                        else{
     349                                mat=NULL;
     350                        }
    297351                        break;
    298352
    … …  
    350404                Vec* pvec=(Vec*)pvalue;
    351405                Vec  outvec=NULL;
    352                 VecDuplicate(vec,&outvec);
    353                 VecCopy(vec,outvec);
     406                if(vec){
     407                        VecDuplicate(vec,&outvec);
     408                        VecCopy(vec,outvec);
     409                }
    354410                *pvec=outvec;
    355411        }
    … …  
    357413                Mat* pmat=(Mat*)pvalue;
    358414                Mat  outmat=NULL;
    359                 MatDuplicate(mat,MAT_COPY_VALUES,&outmat);
     415                if(mat){
     416                        MatDuplicate(mat,MAT_COPY_VALUES,&outmat);
     417                }
    360418                *pmat=outmat;
    361419        }
    … …  
    424482
    425483}
     484
     485#undef __FUNCT__
     486#define __FUNCT__ "SetVec"
     487void  Param::SetVec(Vec value){
     488
     489        if(value){
     490                VecDuplicate(value,&vec);
     491                VecCopy(value,vec);
     492                VecGetSize(value,&M);
     493                N=1;
     494        }
     495        else{
     496                vec=NULL;
     497                M=0;
     498                N=0;
     499        }
     500
     501}

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

    r1 r202  
    4646                void  SetDouble(double value);
    4747                void  SetDoubleVec(double* value,int size);
     48                void  SetVec(Vec value);
    4849                void  SetInteger(int value);
    4950                void  SetString(char* value);

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

    r177 r202  
    663663                GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
    664664
    665                 DL_scalar=alpha2*gauss_weight*Jdet;
     665                if (velocity_param){
     666                        DL_scalar=alpha2*gauss_weight*Jdet;
     667                }
     668                else{
     669                        DL_scalar=0;
     670                }
     671                       
    666672                for (i=0;i<2;i++){
    667673                        DL[i][i]=DL_scalar;
    … …  
    916922        B_param=ParameterInputsRecover(inputs,"B");
    917923        if(temperature_average_param){
    918                 temperature_average-0;
    919924                for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
    920925                temperature_average=temperature_average/3.0;
    … …  
    21222127                                        &Ke_gg_gaussian[0][0],0);
    21232128
    2124                
    21252129                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    21262130                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    … …  
    23522356                /*Get bed slope: */
    23532357                GetParameterDerivativeValue(&slope[0], &b[0],&xyz_list[0][0], gauss_l1l2l3);
    2354                 dbdx=-slope[0];
    2355                 dbdy=-slope[1];
     2358                dbdx=slope[0];
     2359                dbdy=slope[1];
    23562360
    23572361                /* Get Jacobian determinant: */

  • issm/trunk/src/m/classes/public/setelementstype.m

    r106 r202  
    119119deadgrids(find(md.gridonbed))=0;                      %grids from elements on bed are not dead
    120120md.deadgrids=deadgrids;
     121
     122%figure out solution types
     123md.ishutter=any(md.elements_type(:,1)==hutterenum());
     124md.ismacayealpattyn=any(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum() );
     125md.isstokes=any(md.elements_type(:,2));
     126
    121127end

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

    r105 r202  
    1313        disp('   generating degrees of freedom...');
    1414        [m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,parameters);
    15        
     15
    1616        disp('   generating single point constraints...');
    1717        [m.nodes,m.yg]=SpcNodes(m.nodes,m.constraints);
    … …  
    3131        disp('   configuring element and loads...');
    3232        [m.elements,m.loads,m.nodes] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials);
    33        
     33
    3434        disp('   processing parameters...');
    3535        parameters= ProcessParams(parameters,m.part);

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

    r163 r202  
    1818        md.dof=m_dh.nodesets.fsize; %biggest dof number
    1919
     20        %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
     28
    2029        %Get horizontal solution.
    2130        disp(sprintf('\n%s',['computing horizontal velocities...']));
    22 
    23         %plug existing velocity in inputs
    24         if ~isnan(md.vx) & ~isnan(md.vy),
    25                 indx=m_dh.part(1:2:end);
    26                 indy=m_dh.part(2:2:end);
    27                 velocity=zeros(m_dh.nodesets.gsize,1);
    28                 velocity(indx)=md.vx/md.yts;
    29                 velocity(indy)=md.vy/md.yts;
    30                 inputs=struct('velocity',velocity);
    31         else
    32                 inputs={};
    33         end
    3431
    3532        u_g=diagnostic_core_nonlinear(m_dh,inputs);

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

    r133 r202  
    1111        count=1;
    1212        converged=0;
     13
    1314        [velocity_param velocity_is_present]=recover_input(inputs,'velocity');
    1415        if velocity_is_present
    … …  
    1819        end
    1920        soln(count).u_f=[];
    20 
    2121
    2222        if m.parameters.debug,
    … …  
    7474        %   Figure out if convergence is reached.
    7575                if(count>=3 | velocity_is_present),
    76                         dug=soln(count).u_g-soln(count-1).u_g; nduinf=norm(dug,inf)*m.parameters.yts; ndu=norm(dug,2); nu=norm(soln(count-1).u_g,2);
     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);
    7780
    7881                        %Relative criterion

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

    r34 r202  
    1717%Create fem structure (input of diagnostic3d)
    1818fem=struct();
     19
    1920%Figure out which type of elements are present
    20 [fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type);
     21fem.ishutter=ishutter;
     22fem.ismacayealpattyn=ismacayealpattyn;
     23fem.isstokes=isstokes;
    2124
    2225if strcmpi(md.type,'2d'),

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

    r117 r202  
    3636        /*Create constraints: */
    3737        CreateConstraints(&constraints,model,MODEL);
    38        
     38
    3939        /*Create loads: */
    4040        CreateLoads(&loads,model,MODEL);
    41        
     41
    4242        /*Create parameters: */
     43
    4344        CreateParameters(&parameters,model,MODEL);
    4445
Note: See TracChangeset for help on using the changeset viewer.