Changeset 1046


Ignore:
Timestamp:
06/22/09 14:00:59 (15 years ago)
Author:
Mathieu Morlighem
Message:

param_g is now on 1 dof

Location:
issm/trunk/src
Files:
15 edited

Legend:

Unmodified
Added
Removed

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

    r1 r1046  
    1313#include "../EnumDefinitions/EnumDefinitions.h"
    1414
    15 void ControlConstrainx( double* p_g, int gsize, double mincontrolconstraint, double maxcontrolconstraint,char* control_type){
     15void ControlConstrainx( double* p_g, int numdofnodes, double mincontrolconstraint, double maxcontrolconstraint,char* control_type){
    1616
    1717        int i;
    … …  
    2121        }
    2222        else{
    23                 for(i=0;i<gsize;i=i+2){
     23                for(i=0;i<numdofnodes;i++){
    2424
    2525                        if(!isnan(mincontrolconstraint)){

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

    r1 r1046  
    99
    1010/* local prototypes: */
    11 void ControlConstrainx( double* p_g, int gsize, double mincontrolconstraint, double maxcontrolconstraint,char* control_type);
     11void ControlConstrainx( double* p_g, int numberofnodes, double mincontrolconstraint, double maxcontrolconstraint,char* control_type);
    1212
    1313#endif  /* _CONTROLCONSTRAINX_H */

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

    r465 r1046  
    1313#include "../EnumDefinitions/EnumDefinitions.h"
    1414
    15 void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
     15void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    1616                double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    17 
    18         int i;
    19         int gsize;
    20         int found;
    2117
    2218        /*output: */
    … …  
    2622        elements->Configure(elements,loads, nodes, materials);
    2723
    28         /*Get size of matrix: */
    29         gsize=nodes->NumberOfDofs();
    30 
    3124        /*Allocate grad_g: */
    32         grad_g=NewVec(gsize);
     25        grad_g=NewVec(numberofnodes);
    3326
    3427        /*Compute gradients: */

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

    r465 r1046  
    99
    1010/* local prototypes: */
    11 void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
     11void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,
    1212                double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type);
    1313

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

    r774 r1046  
    135135        parameters->AddObject(param);
    136136       
    137         param_g=(double*)xcalloc(model->numberofnodes*2,sizeof(double));
    138         for(i=0;i<model->numberofnodes;i++)param_g[2*i+0]=control_parameter[i];
     137        param_g=(double*)xcalloc(model->numberofnodes,sizeof(double));
     138        for(i=0;i<model->numberofnodes;i++)param_g[i]=control_parameter[i];
    139139
    140140        count++;
    141141        param= new Param(count,"param_g",DOUBLEVEC);
    142         param->SetDoubleVec(param_g,2*model->numberofnodes,2);
     142        param->SetDoubleVec(param_g,model->numberofnodes,1);
    143143        parameters->AddObject(param);
    144144

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

    r983 r1046  
    44
    55#ifdef HAVE_CONFIG_H
    6         #include "config.h"
     6#include "config.h"
    77#else
    88#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    … …  
    3333
    3434Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3],
    35                                 double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
    36                                 double tria_viscosity_overshoot,int tria_artdiff){
    37        
     35                        double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
     36                        double tria_viscosity_overshoot,int tria_artdiff){
     37
    3838        int i;
    39        
     39
    4040        id=tria_id;
    4141        mid=tria_mid;
    … …  
    6666        viscosity_overshoot=tria_viscosity_overshoot;
    6767        artdiff=tria_artdiff;
    68        
     68
    6969        return;
    7070}
    … …  
    156156        /*get enum type of Tria: */
    157157        enum_type=TriaEnum();
    158        
     158
    159159        /*marshall enum: */
    160160        memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
    161        
     161
    162162        /*marshall Tria data: */
    163163        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
    … …  
    187187        memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
    188188        memcpy(marshalled_dataset,&artdiff,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff);
    189        
     189
    190190        *pmarshalled_dataset=marshalled_dataset;
    191191        return;
    192192}
    193                
     193
    194194int   Tria::MarshallSize(){
    195195        return sizeof(id)
    196                 +sizeof(mid)
    197                 +sizeof(mparid)
    198                 +sizeof(node_ids)
    199                 +sizeof(nodes)
    200                 +sizeof(node_offsets)
    201                 +sizeof(matice)
    202                 +sizeof(matice_offset)
    203                 +sizeof(matpar)
    204                 +sizeof(matpar_offset)
    205                 +sizeof(h)
    206                 +sizeof(s)
    207                 +sizeof(b)
    208                 +sizeof(k)
    209                 +sizeof(melting)
    210                 +sizeof(accumulation)
    211                 +sizeof(geothermalflux)
    212                 +sizeof(friction_type)
    213                 +sizeof(onbed)
    214                 +sizeof(p)
    215                 +sizeof(q)
    216                 +sizeof(shelf)
    217                 +sizeof(meanvel)
    218                 +sizeof(epsvel)
    219                 +sizeof(viscosity_overshoot)
    220                 +sizeof(artdiff)
    221                 +sizeof(int); //sizeof(int) for enum type
     196          +sizeof(mid)
     197          +sizeof(mparid)
     198          +sizeof(node_ids)
     199          +sizeof(nodes)
     200          +sizeof(node_offsets)
     201          +sizeof(matice)
     202          +sizeof(matice_offset)
     203          +sizeof(matpar)
     204          +sizeof(matpar_offset)
     205          +sizeof(h)
     206          +sizeof(s)
     207          +sizeof(b)
     208          +sizeof(k)
     209          +sizeof(melting)
     210          +sizeof(accumulation)
     211          +sizeof(geothermalflux)
     212          +sizeof(friction_type)
     213          +sizeof(onbed)
     214          +sizeof(p)
     215          +sizeof(q)
     216          +sizeof(shelf)
     217          +sizeof(meanvel)
     218          +sizeof(epsvel)
     219          +sizeof(viscosity_overshoot)
     220          +sizeof(artdiff)
     221          +sizeof(int); //sizeof(int) for enum type
    222222}
    223223
    … …  
    291291
    292292        int i;
    293        
     293
    294294        DataSet* loadsin=NULL;
    295295        DataSet* nodesin=NULL;
    … …  
    303303        /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
    304304        ResolvePointers((Object**)nodes,node_ids,node_offsets,3,nodesin);
    305        
     305
    306306        /*Same for materials: */
    307307        ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
    … …  
    317317        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    318318        if (analysis_type==ControlAnalysisEnum()){
    319                
     319
    320320                CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
    321321        }
    322322        else if (analysis_type==DiagnosticAnalysisEnum()){
    323        
     323
    324324                if (sub_analysis_type==HorizAnalysisEnum()){
    325325
    … …  
    361361        int          doflist[numdof];
    362362        int          numberofdofspernode;
    363        
     363
    364364        /* gaussian points: */
    365365        int     num_gauss,ig;
    … …  
    375375        double newviscosity; //viscosity
    376376        double oldviscosity; //viscosity
    377        
     377
    378378        /* strain rate: */
    379379        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    … …  
    389389        double Ke_gg[numdof][numdof]; //local element stiffness matrix
    390390        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    391        
     391
    392392        double Jdet;
    393        
     393
    394394        /*input parameters for structural analysis (diagnostic): */
    395395        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    … …  
    414414        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    415415
    416         #ifdef _DEBUGELEMENTS_
     416#ifdef _DEBUGELEMENTS_
    417417        if(my_rank==RANK && id==ELID){
    418418                printf("El id %i Rank %i TriaElemnet input list before gaussian loop: \n",ELID,RANK);
    419419                printf("   rho_ice: %g \n",matpar->GetRhoIce());
    420420                printf("   gravity: %g \n",matpar->GetG())
    421                 printf("   rho_water: %g \n",matpar->GetRhoWater());
     421                  printf("   rho_water: %g \n",matpar->GetRhoWater());
    422422                printf("   Velocity: \n");
    423423                for (i=0;i<numgrids;i++){
    … …  
    430430                printf("   bed [%g %g %g]\n",b[0],b[1],b[2]);
    431431        }
    432         #endif
     432#endif
    433433
    434434
    … …  
    436436        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    437437
    438         #ifdef _DEBUGELEMENTS_
     438#ifdef _DEBUGELEMENTS_
    439439        if(my_rank==RANK && id==ELID){
    440440                printf("   gaussian points: \n");
    … …  
    443443                }
    444444        }
    445         #endif
     445#endif
    446446
    447447        /* Start  looping on the number of gaussian points: */
    … …  
    460460                GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
    461461                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
    462                
     462
    463463                /*Get viscosity: */
    464464                matice->GetViscosity2d(&viscosity, &epsilon[0]);
    465465                matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
    466                
     466
    467467                /* Get Jacobian determinant: */
    468468                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    469469
    470470                /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    471                    onto this scalar matrix, so that we win some computational time: */
     471                        onto this scalar matrix, so that we win some computational time: */
    472472                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    473473                D_scalar=newviscosity*thickness*gauss_weight*Jdet;
    … …  
    476476                        D[i][i]=D_scalar;
    477477                }
    478                
    479                 #ifdef _DEBUGELEMENTS_
     478
     479#ifdef _DEBUGELEMENTS_
    480480                if(my_rank==RANK && id==ELID){
    481481                        printf("   gaussian loop %i\n",ig);
    … …  
    492492                        printf("      gauss_weight: %g \n",gauss_weight);
    493493                }
    494                 #endif
     494#endif
    495495
    496496                /*Get B and Bprime matrices: */
    … …  
    500500                /*  Do the triple product tB*D*Bprime: */
    501501                TripleMultiply( &B[0][0],3,numdof,1,
    502                                           &D[0][0],3,3,0,
    503                                           &Bprime[0][0],3,numdof,0,
    504                                           &Ke_gg_gaussian[0][0],0);
     502                                        &D[0][0],3,3,0,
     503                                        &Bprime[0][0],3,numdof,0,
     504                                        &Ke_gg_gaussian[0][0],0);
    505505
    506506                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    507507                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    508                
    509                 #ifdef _DEBUGELEMENTS_
     508
     509#ifdef _DEBUGELEMENTS_
    510510                if(my_rank==RANK && id==ELID){
    511511                        printf("      B:\n");
    … …  
    524524                        }
    525525                }
    526                 #endif
     526#endif
    527527        } // for (ig=0; ig<num_gauss; ig++)
    528528
    … …  
    537537
    538538
    539         #ifdef _DEBUGELEMENTS_
     539#ifdef _DEBUGELEMENTS_
    540540        if(my_rank==RANK && id==ELID){
    541541                printf("      Ke_gg erms:\n");
    … …  
    552552
    553553        }
    554         #endif
    555 
    556         cleanup_and_return:
     554#endif
     555
     556cleanup_and_return:
    557557        xfree((void**)&first_gauss_area_coord);
    558558        xfree((void**)&second_gauss_area_coord);
    … …  
    576576        int          doflist[numdof];
    577577        int          numberofdofspernode;
    578        
     578
    579579        /* gaussian points: */
    580580        int     num_gauss,ig;
    … …  
    597597        double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    598598        double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    599        
     599
    600600        double Jdettria;
    601        
     601
    602602        /*input parameters for structural analysis (diagnostic): */
    603603        double  vxvy_list[numgrids][2]={0.0};
    … …  
    628628                vy_list[i]=vxvy_list[i][1];
    629629        }
    630        
     630
    631631        found=inputs->Recover("dt",&dt);
    632632        if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
    … …  
    645645                v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
    646646                v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
    647                
     647
    648648                K[0][0]=pow(Jdettria,.5)/2.0*fabs(v_gauss[0]);
    649649                K[1][1]=pow(Jdettria,.5)/2.0*fabs(v_gauss[1]);
    … …  
    671671                /*  Do the triple product tL*D*L: */
    672672                TripleMultiply( &L[0],1,numdof,1,
    673                                           &DL_scalar,1,1,0,
    674                                           &L[0],1,numdof,0,
    675                                           &Ke_gg_gaussian[0][0],0);
    676                
     673                                        &DL_scalar,1,1,0,
     674                                        &L[0],1,numdof,0,
     675                                        &Ke_gg_gaussian[0][0],0);
     676
    677677                /*Get B  and B prime matrix: */
    678678                GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
    679679                GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
    680                
     680
    681681                //Get vx, vy and their derivatives at gauss point
    682682                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
    683683                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
    684                
     684
    685685                GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
    686686                GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
    … …  
    702702
    703703                TripleMultiply( &B[0][0],2,numdof,1,
    704                                           &DL[0][0],2,2,0,
    705                                           &B[0][0],2,numdof,0,
    706                                           &Ke_gg_thickness1[0][0],0);
     704                                        &DL[0][0],2,2,0,
     705                                        &B[0][0],2,numdof,0,
     706                                        &Ke_gg_thickness1[0][0],0);
    707707
    708708                TripleMultiply( &B[0][0],2,numdof,1,
    709                                           &DLprime[0][0],2,2,0,
    710                                           &Bprime[0][0],2,numdof,0,
    711                                           &Ke_gg_thickness2[0][0],0);
     709                                        &DLprime[0][0],2,2,0,
     710                                        &Bprime[0][0],2,numdof,0,
     711                                        &Ke_gg_thickness2[0][0],0);
    712712
    713713                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    … …  
    715715                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
    716716                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
    717                
     717
    718718                if(artdiff){
    719                        
     719
    720720                        /* Compute artificial diffusivity */
    721721                        KDL[0][0]=DL_scalar*K[0][0];
    … …  
    723723
    724724                        TripleMultiply( &Bprime[0][0],2,numdof,1,
    725                                                   &KDL[0][0],2,2,0,
    726                                                   &Bprime[0][0],2,numdof,0,
    727                                                   &Ke_gg_gaussian[0][0],0);
     725                                                &KDL[0][0],2,2,0,
     726                                                &Bprime[0][0],2,numdof,0,
     727                                                &Ke_gg_gaussian[0][0],0);
    728728
    729729                        /* Add artificial diffusivity matrix */
    730730                        for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    731                        
    732                 }
    733 
    734                 #ifdef _DEBUGELEMENTS_
     731
     732                }
     733
     734#ifdef _DEBUGELEMENTS_
    735735                if(my_rank==RANK && id==ELID){
    736736                        printf("      B:\n");
    … …  
    749749                        }
    750750                }
    751                 #endif
     751#endif
    752752        } // for (ig=0; ig<num_gauss; ig++)
    753753
    … …  
    755755        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    756756
    757         #ifdef _DEBUGELEMENTS_
     757#ifdef _DEBUGELEMENTS_
    758758        if(my_rank==RANK && id==ELID){
    759759                printf("      Ke_gg erms:\n");
    … …  
    770770
    771771        }
    772         #endif
    773 
    774         cleanup_and_return:
     772#endif
     773
     774cleanup_and_return:
    775775        xfree((void**)&first_gauss_area_coord);
    776776        xfree((void**)&second_gauss_area_coord);
    … …  
    795795        int          doflist[numdof];
    796796        int          numberofdofspernode;
    797        
     797
    798798        /* gaussian points: */
    799799        int     num_gauss,ig;
    … …  
    809809        double L[numgrids];
    810810        double Jdettria;
    811        
     811
    812812        /*input parameters for structural analysis (diagnostic): */
    813813        double  accumulation_list[numgrids]={0.0};
    … …  
    861861                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
    862862                GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3);
    863                
     863
    864864                /* Add value into pe_g: */
    865865                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
    866                
     866
    867867        } // for (ig=0; ig<num_gauss; ig++)
    868868
    … …  
    870870        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    871871
    872         cleanup_and_return:
     872cleanup_and_return:
    873873        xfree((void**)&first_gauss_area_coord);
    874874        xfree((void**)&second_gauss_area_coord);
    … …  
    893893        int          doflist[numdof];
    894894        int          numberofdofspernode;
    895        
     895
    896896        /* gaussian points: */
    897897        int     num_gauss,ig;
    … …  
    911911        double Ke_gg[numdof][numdof]; //local element stiffness matrix
    912912        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    913        
     913
    914914        double Jdet;
    915        
     915
    916916        /*slope: */
    917917        double  slope[2]={0.0,0.0};
    … …  
    933933        /*recover pointers: */
    934934        inputs=(ParameterInputs*)vinputs;
    935        
     935
    936936        /* Get node coordinates and dof list: */
    937937        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    953953        /*Build alpha2_list used by drag stiffness matrix*/
    954954        Friction* friction=NewFriction();
    955        
     955
    956956        /*Initialize all fields: */
    957957        friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    958958        strcpy(friction->element_type,"2d");
    959        
     959
    960960        friction->gravity=matpar->GetG();
    961961        friction->rho_ice=matpar->GetRhoIce();
    … …  
    974974        DeleteFriction(&friction);
    975975
    976         #ifdef _DEBUGELEMENTS_
     976#ifdef _DEBUGELEMENTS_
    977977        if(my_rank==RANK && id==ELID){
    978978                printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
    979979        }
    980         #endif
     980#endif
    981981
    982982        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    983983        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    984984
    985         #ifdef _DEBUGELEMENTS_
     985#ifdef _DEBUGELEMENTS_
    986986        if(my_rank==RANK && id==ELID){
    987987                printf("   gaussian points: \n");
    … …  
    990990                }
    991991        }
    992         #endif
     992#endif
    993993
    994994        /* Start  looping on the number of gaussian points: */
    … …  
    10251025                        DL[i][i]=DL_scalar;
    10261026                }
    1027                
     1027
    10281028                /*  Do the triple producte tL*D*L: */
    10291029                TripleMultiply( &L[0][0],2,numdof,1,
    … …  
    10391039        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    10401040
    1041         cleanup_and_return:
     1041cleanup_and_return:
    10421042        xfree((void**)&first_gauss_area_coord);
    10431043        xfree((void**)&second_gauss_area_coord);
    … …  
    10621062        int          doflist[numdof];
    10631063        int          numberofdofspernode;
    1064        
     1064
    10651065        /* gaussian points: */
    10661066        int     num_gauss,ig;
    … …  
    10791079        double Ke_gg[numdof][numdof]; //local element stiffness matrix
    10801080        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    1081        
     1081
    10821082        double Jdet;
    1083        
     1083
    10841084        /* Get node coordinates and dof list: */
    10851085        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    11001100                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    11011101
    1102                
     1102
    11031103                /*Get L matrix: */
    11041104                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
    … …  
    11061106                /* Get Jacobian determinant: */
    11071107                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1108                
     1108
    11091109                DL_scalar=gauss_weight*Jdet;
    11101110
    … …  
    11211121        /*Add Ke_gg to global matrix Kgg: */
    11221122        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    1123                
    1124         cleanup_and_return:
     1123
     1124cleanup_and_return:
    11251125        xfree((void**)&first_gauss_area_coord);
    11261126        xfree((void**)&second_gauss_area_coord);
    … …  
    11321132#define __FUNCT__ "Tria::CreatePVector"
    11331133void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
    1134        
     1134
    11351135        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    11361136        if (analysis_type==ControlAnalysisEnum()){
    1137                
     1137
    11381138                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
    1139        
     1139
    11401140        }
    11411141        else if (analysis_type==DiagnosticAnalysisEnum()){
    11421142                if (sub_analysis_type==HorizAnalysisEnum()){
    1143                
     1143
    11441144                        CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
    1145                
     1145
    11461146                }
    11471147                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
    11481148        }
    11491149        else if (analysis_type==SlopeComputeAnalysisEnum()){
    1150                
     1150
    11511151                CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
    11521152        }
    … …  
    11741174        int          doflist[numdof];
    11751175        int          numberofdofspernode;
    1176        
     1176
    11771177        /* parameters: */
    11781178        double  plastic_stress;
    … …  
    12151215
    12161216
    1217         #ifdef _DEBUGELEMENTS_
     1217#ifdef _DEBUGELEMENTS_
    12181218        if(my_rank==RANK && id==ELID){
    12191219                printf("gravity %g\n",matpar->GetG());
    … …  
    12241224                printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]);
    12251225        }
    1226         #endif
     1226#endif
    12271227
    12281228
    … …  
    12301230        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
    12311231
    1232         #ifdef _DEBUGELEMENTS_
     1232#ifdef _DEBUGELEMENTS_
    12331233        if(my_rank==RANK && id==ELID){
    12341234                printf("   gaussian points: \n");
    … …  
    12371237                }
    12381238        }
    1239         #endif
     1239#endif
    12401240
    12411241
    … …  
    12511251                /*Compute thickness at gaussian point: */
    12521252                GetParameterValue(&thickness, &h[0],gauss_l1l2l3);
    1253        
     1253
    12541254                GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
    1255                
     1255
    12561256                /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the
    12571257                 * element itself: */
    … …  
    12621262                /* Get Jacobian determinant: */
    12631263                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1264                
    1265                  /*Get nodal functions: */
     1264
     1265                /*Get nodal functions: */
    12661266                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    12671267
    … …  
    12701270
    12711271
    1272                 #ifdef _DEBUGELEMENTS_
     1272#ifdef _DEBUGELEMENTS_
    12731273                if(my_rank==RANK && id==ELID){
    12741274                        printf("      gaussian %i\n",ig);
    … …  
    12801280                        if(friction_type==1)printf("      plastic_stress(%g)\n",plastic_stress);
    12811281                }
    1282                 #endif
     1282#endif
    12831283
    12841284                /*Build pe_g_gaussian vector: */
    … …  
    13031303        } //for (ig=0; ig<num_gauss; ig++)
    13041304
    1305         #ifdef _DEBUGELEMENTS_
     1305#ifdef _DEBUGELEMENTS_
    13061306        if(my_rank==RANK && id==ELID){
    13071307                printf("      pe_g->terms\n",ig);
    … …  
    13151315                }
    13161316        }
    1317         #endif
     1317#endif
    13181318
    13191319        /*Add pe_g to global vector pg: */
    13201320        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    13211321
    1322         cleanup_and_return:
     1322cleanup_and_return:
    13231323        xfree((void**)&first_gauss_area_coord);
    13241324        xfree((void**)&second_gauss_area_coord);
    … …  
    13421342        int          doflist[numdof];
    13431343        int          numberofdofspernode;
    1344        
     1344
    13451345        /* gaussian points: */
    13461346        int     num_gauss,ig;
    … …  
    13911391
    13921392                GetParameterDerivativeValue(&slope[0], &param[0],&xyz_list[0][0], gauss_l1l2l3);
    1393                
     1393
    13941394                /* Get Jacobian determinant: */
    13951395                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1396                
    1397                  /*Get nodal functions: */
     1396
     1397                /*Get nodal functions: */
    13981398                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    13991399
    … …  
    14141414        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    14151415
    1416         cleanup_and_return:
     1416cleanup_and_return:
    14171417        xfree((void**)&first_gauss_area_coord);
    14181418        xfree((void**)&second_gauss_area_coord);
    … …  
    14451445        inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes);
    14461446        inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes);
    1447        
     1447
    14481448        //Update material if necessary
    14491449        if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){
    … …  
    14521452                matice->SetB(B_average);
    14531453        }
    1454        
     1454
    14551455        if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){
    14561456                B_average=(B_list[0]+B_list[1]+B_list[2])/3.0;
    … …  
    14591459
    14601460}
    1461                
     1461
    14621462void  Tria::GetDofList(int* doflist,int* pnumberofdofspernode){
    14631463
    … …  
    14651465        int doflist_per_node[MAXDOFSPERNODE];
    14661466        int numberofdofspernode;
    1467        
     1467
    14681468        for(i=0;i<3;i++){
    14691469                nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
    … …  
    14911491#define __FUNCT__ "Tria::GetParameterValue"
    14921492void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){
    1493        
     1493
    14941494        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian
    14951495         * point specifie by gauss_l1l2l3: */
    1496        
     1496
    14971497        /*nodal functions: */
    14981498        double l1l2l3[3];
    … …  
    15131513#define __FUNCT__ "Tria::GetParameterDerivativeValue"
    15141514void Tria::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3){
    1515          
     1515
    15161516        const int NDOF2=2;
    15171517        const int numgrids=3;
    … …  
    15231523         * p is a vector of size 2x1 already allocated.
    15241524         */
    1525        
     1525
    15261526        double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system.
    15271527
    … …  
    15481548        GetB(&B[0][0], xyz_list, gauss_l1l2l3);
    15491549
    1550         #ifdef _DEBUG_
     1550#ifdef _DEBUG_
    15511551        printf("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
    15521552        printf("              [ %lf   %lf  \n",B[1][0],B[1][1]);
    … …  
    15601560        printf("              [ %lf   %lf  \n",B[1][4],B[1][5]);
    15611561        printf("              [ %lf   %lf  ]\n",B[2][4],B[2][5]);
    1562                
     1562
    15631563        for (i=0;i<numgrids;i++){
    15641564                printf("Velocity for grid %i %lf %lf\n",i,*(vxvy_list+2*i+0),*(vxvy_list+2*i+1));
    15651565        }
    1566         #endif
     1566#endif
    15671567
    15681568        /*Multiply B by velocity, to get strain rate: */
    15691569        MatrixMultiply( &B[0][0],3,NDOF2*numgrids,0,
    1570                                       velocity,NDOF2*numgrids,1,0,
    1571                                                   epsilon,0);
     1570                                velocity,NDOF2*numgrids,1,0,
     1571                                epsilon,0);
    15721572
    15731573}
    … …  
    15811581
    15821582        double x1,x2,x3,y1,y2,y3;
    1583        
     1583
    15841584        x1=*(xyz_list+3*0+0);
    15851585        y1=*(xyz_list+3*0+1);
    … …  
    15961596                printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!");
    15971597        }
    1598        
     1598
    15991599}
    16001600
    … …  
    16071607
    16081608        double x1,x2,x3,y1,y2,y3,z1,z2,z3;
    1609        
     1609
    16101610        x1=*(xyz_list+3*0+0);
    16111611        y1=*(xyz_list+3*0+1);
    … …  
    16251625                printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!");
    16261626        }
    1627        
     1627
    16281628}
    16291629
    … …  
    16431643         * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
    16441644         */
    1645        
     1645
    16461646        int i;
    16471647        const int NDOF2=2;
    … …  
    16541654        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);
    16551655
    1656         #ifdef _DEBUG_
     1656#ifdef _DEBUG_
    16571657        for (i=0;i<3;i++){
    16581658                printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh2dh3_basic[0][i],dh1dh2dh3_basic[1][i]);
    16591659        }
    1660         #endif
     1660#endif
    16611661
    16621662        /*Build B: */
    … …  
    16861686         * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
    16871687         */
    1688        
     1688
    16891689        int i;
    16901690        const int NDOF2=2;
    … …  
    17261726         * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
    17271727         */
    1728        
     1728
    17291729        int i;
    17301730        const int NDOF2=2;
    … …  
    17371737        GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    17381738
    1739         #ifdef _DELUG_
     1739#ifdef _DELUG_
    17401740        for (i=0;i<3;i++){
    17411741                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
    17421742        }
    1743         #endif
     1743#endif
    17441744
    17451745        /*Build L: */
    … …  
    17731773         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
    17741774         */
    1775        
     1775
    17761776        int i;
    17771777        const int NDOF1=1;
    … …  
    17841784        GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
    17851785
    1786         #ifdef _DEBUG_
     1786#ifdef _DEBUG_
    17871787        for (i=0;i<3;i++){
    17881788                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
    17891789        }
    1790         #endif
     1790#endif
    17911791
    17921792        /*Build B_prog: */
    … …  
    18111811         * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
    18121812         */
    1813        
     1813
    18141814        int i;
    18151815        const int NDOF1=1;
    … …  
    18341834#define __FUNCT__ "Tria::GetNodalFunctions"
    18351835void Tria::GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3){
    1836        
     1836
    18371837        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    18381838
    … …  
    18511851#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesBasic"
    18521852void Tria::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3){
    1853        
     1853
    18541854        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    18551855         * basic coordinate system: */
    … …  
    18851885#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesParams"
    18861886void Tria::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3){
    1887        
     1887
    18881888        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    18891889         * natural coordinate system) at the gaussian point. */
    … …  
    19151915        const int NDOF2=2;
    19161916        const int numgrids=3;
    1917        
     1917
    19181918        /*Call Jacobian routine to get the jacobian:*/
    19191919        GetJacobian(Jinv, xyz_list, gauss_l1l2l3);
    … …  
    19211921        /*Invert Jacobian matrix: */
    19221922        MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
    1923                
     1923
    19241924}
    19251925
    … …  
    19351935        double x1,y1,x2,y2,x3,y3;
    19361936        double sqrt3=sqrt(3.0);
    1937        
     1937
    19381938        x1=*(xyz_list+numgrids*0+0);
    19391939        y1=*(xyz_list+numgrids*0+1);
    … …  
    19581958}
    19591959
    1960                
     1960
    19611961void  Tria::GetNodes(void** vpnodes){
    19621962        int i;
    … …  
    19671967        }
    19681968}
    1969                
     1969
    19701970
    19711971int Tria::GetOnBed(){
    … …  
    19791979}
    19801980void          Tria::GetBedList(double* bed_list){
    1981        
     1981
    19821982        int i;
    19831983        for(i=0;i<3;i++)bed_list[i]=b[i];
    19841984
    19851985}
    1986  
     1986
    19871987
    19881988Object* Tria::copy() {
    1989        
     1989
    19901990        return new Tria(*this);
    19911991
    … …  
    19981998
    19991999        int i;
    2000        
     2000
    20012001        /* node data: */
    20022002        const int    numgrids=3;
    … …  
    21042104                throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
    21052105        }
    2106        
     2106
    21072107        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    21082108        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    21092109
    2110         #ifdef _DEBUGELEMENTS_
     2110#ifdef _DEBUGELEMENTS_
    21112111        if(my_rank==RANK && id==ELID){
    21122112                printf("   gaussian points: \n");
    … …  
    21152115                }
    21162116        }
    2117         #endif
    2118        
     2117#endif
     2118
    21192119        /* Start  looping on the number of gaussian points: */
    21202120        for (ig=0; ig<num_gauss; ig++){
    … …  
    21272127                /* Get Jacobian determinant: */
    21282128                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    2129                 #ifdef _DEBUG_
     2129#ifdef _DEBUG_
    21302130                printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
    2131                 #endif
    2132                
     2131#endif
     2132
    21332133                /* Get nodal functions value at gaussian point:*/
    21342134                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    … …  
    21782178                        throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
    21792179                }
    2180                
     2180
    21812181                /*Add due_g_gaussian vector to due_g: */
    21822182                for( i=0; i<numdof; i++){
    … …  
    21842184                }
    21852185        }
    2186        
     2186
    21872187        /*Add due_g to global vector du_g: */
    21882188        VecSetValues(du_g,numdof,doflist,(const double*)due_g,ADD_VALUES);
    2189        
    2190         cleanup_and_return:
     2189
     2190cleanup_and_return:
    21912191        xfree((void**)&first_gauss_area_coord);
    21922192        xfree((void**)&second_gauss_area_coord);
    … …  
    22212221        double       xyz_list[numgrids][3];
    22222222        int          doflist[numdof];
     2223        int          doflist1[numgrids];
    22232224        int          numberofdofspernode;
    22242225
    … …  
    22512252        double  lambda,mu;
    22522253        double  bed,thickness,Neff;
    2253        
     2254
    22542255        /*drag: */
    22552256        double  pcoeff,qcoeff;
    … …  
    22592260
    22602261        /*element vector at the gaussian points: */
    2261         double  grade_g[numdof];
     2262        double  grade_g[numgrids];
    22622263        double  grade_g_gaussian[numgrids];
    22632264
    … …  
    22792280        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    22802281        GetDofList(&doflist[0],&numberofdofspernode);
     2282        GetDofList1(&doflist1[0]);
    22812283
    22822284        /* Set grade_g to 0: */
    2283         for(i=0;i<numdof;i++) grade_g[i]=0.0;
     2285        for(i=0;i<numgrids;i++) grade_g[i]=0.0;
    22842286
    22852287        /* recover input parameters: */
    … …  
    23882390               
    23892391                /*Add gradje_g_gaussian vector to gradje_g: */
    2390                 for( i=0; i<numgrids; i++)grade_g[i*numberofdofspernode]+=grade_g_gaussian[i];
     2392                for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
    23912393        }
    23922394
    23932395        /*Add grade_g to global vector grad_g: */
    2394         VecSetValues(grad_g,numdof,doflist,(const double*)grade_g,ADD_VALUES);
     2396        VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    23952397       
    23962398        cleanup_and_return:
    … …  
    24152417        double       xyz_list[numgrids][3];
    24162418        int          doflist[numdof];
     2419        int          doflist1[numgrids];
    24172420        int          numberofdofspernode;
    24182421
    … …  
    24342437
    24352438        /*element vector at the gaussian points: */
    2436         double  grade_g[numdof];
     2439        double  grade_g[numgrids];
    24372440        double  grade_g_gaussian[numgrids];
    24382441
    … …  
    24662469        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    24672470        GetDofList(&doflist[0],&numberofdofspernode);
     2471        GetDofList1(&doflist1[0]);
    24682472
    24692473        /* Set grade_g to 0: */
    2470         for(i=0;i<numdof;i++) grade_g[i]=0.0;
     2474        for(i=0;i<numgrids;i++) grade_g[i]=0.0;
    24712475
    24722476        /* recover input parameters: */
    … …  
    25292533
    25302534                /*Add grade_g_gaussian to grade_g: */
    2531                 for( i=0; i<numgrids;i++) grade_g[i*numberofdofspernode]+=grade_g_gaussian[i];
     2535                for( i=0; i<numgrids;i++) grade_g[i]+=grade_g_gaussian[i];
    25322536        }
    25332537       
    25342538        /*Add grade_g to global vector grad_g: */
    2535         VecSetValues(grad_g,numdof,doflist,(const double*)grade_g,ADD_VALUES);
     2539        VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    25362540       
    25372541        cleanup_and_return:

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

    r465 r1046  
    2121        int analysis_type;
    2222        int sub_analysis_type;
     23        int numberofnodes;
    2324        char* solverstring=NULL;
    2425        char* control_type=NULL;
    … …  
    4445        femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type");
    4546        femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
     47        femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    4648        femmodel->parameters->FindParam((void*)&solverstring,"solverstring");
    4749        femmodel->parameters->FindParam((void*)&control_type,"control_type");
    … …  
    6971        VecToMPISerial(&lambda_g_double,lambda_g);VecFree(&lambda_g);
    7072
    71         Gradjx( &grad_g, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials,
     73        Gradjx( &grad_g, numberofnodes,femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials,
    7274                u_g_double,lambda_g_double, inputs,analysis_type,sub_analysis_type,control_type);
    7375       

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

    r1037 r1046  
    399399
    400400                        for(i=0;i<numberofnodes;i++){
    401                                 parameter[i]=param_g[2*(int)partition[i]];
     401                                parameter[i]=param_g[(int)partition[i]];
    402402                        }
    403403                       

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

    r1037 r1046  
    9999
    100100                _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
    101                 inputs->Add(control_type,param_g,2,numberofnodes);
     101                inputs->Add(control_type,param_g,1,numberofnodes);
    102102                inputs->Add("fit",fit[n]);
    103103
    … …  
    123123
    124124                _printf_("%s\n","      updating parameter using optimized search scalar...");
    125                 for(i=0;i<gsize;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
     125                for(i=0;i<numberofnodes;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
    126126                _printf_("%s\n","      done.");
    127127
    128128                _printf_("%s\n","      constraining the new distribution...");   
    129                 ControlConstrainx( param_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
     129                ControlConstrainx(param_g,numberofnodes,mincontrolconstraint,maxcontrolconstraint,control_type);
    130130                _printf_("%s\n","      done.");
    131131
    … …  
    138138        /*Write results to disk: */
    139139        _printf_("%s\n","      preparing final velocity solution...");
    140         inputs->Add(control_type,param_g,2,numberofnodes);
     140        inputs->Add(control_type,param_g,1,numberofnodes);
    141141        inputs->Add("fit",fit[n]);
    142142

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

    r758 r1046  
    6464
    6565        /*First copy param_g so we don't modify it: */
    66         param_g_copy=(double*)xmalloc(gsize*sizeof(double));
    67         memcpy(param_g_copy,param_g,gsize*sizeof(double));
     66        param_g_copy=(double*)xmalloc(numberofnodes*sizeof(double));
     67        memcpy(param_g_copy,param_g,numberofnodes*sizeof(double));
    6868
    6969        /*First, update param_g using search_scalar: */
    70         for(i=0;i<gsize;i++)param_g_copy[i]=param_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
     70        for(i=0;i<numberofnodes;i++)param_g_copy[i]=param_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
    7171
    7272        /*Constrain:*/
    73         ControlConstrainx( param_g_copy,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
     73        ControlConstrainx(param_g_copy,numberofnodes,mincontrolconstraint,maxcontrolconstraint,control_type);
    7474
    7575        /*Add new parameter to inputs: */
    76         inputs->Add(control_type,param_g_copy,2,numberofnodes);
     76        inputs->Add(control_type,param_g_copy,1,numberofnodes);
    7777
    7878        //Run diagnostic with updated parameters.
    … …  
    8585                u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type);
    8686
    87 
    8887        /*Free ressources:*/
    8988        xfree((void**)&optscal);

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

    r1033 r1046  
    1818%initialize control parameters, gradients and observations
    1919u_g_obs=m_dh.parameters.u_g_obs;
    20 grad_g=zeros(m_dh.nodesets.gsize,1);
     20grad_g=zeros(m_dh.parameters.numberofnodes,1);
    2121param_g=models.dh.parameters.param_g;
    2222
    … …  
    3333        %update inputs with new fit
    3434        inputs=add(inputs,'fit',m_dh.parameters.fit(n),'double');
    35         inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',2,m_dh.parameters.numberofnodes);
     35        inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',1,m_dh.parameters.numberofnodes);
    3636
    3737        %Update inputs in datasets

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

    r980 r1046  
    66
    77%From parameters, build inputs for icediagnostic_core, using the final parameters
    8 inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes);
     8inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',1,m.parameters.numberofnodes);
    99results.u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type);
    1010

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

    r465 r1046  
    1010
    1111%Plug parameter into inputs
    12 inputs=add(inputs,m.parameters.control_type,parameter,'doublevec',2,m.parameters.numberofnodes);
     12inputs=add(inputs,m.parameters.control_type,parameter,'doublevec',1,m.parameters.numberofnodes);
    1313
    1414%Run diagnostic with updated parameters.

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

    r980 r1046  
    124124                        newresults(i).step=results(i).step;
    125125                        newresults(i).time=results(i).time;
    126                         newresults(i).parameter=param_g(1:2:end);
     126                        newresults(i).parameter=param_g;
    127127
    128128                else

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

    r465 r1046  
    1919        double*  lambda_g=NULL;
    2020        ParameterInputs* inputs=NULL;
    21         char*             analysis_type_string=NULL;
    22         int               analysis_type;
    23         char*             sub_analysis_type_string=NULL;
    24         int               sub_analysis_type;
     21        char*    analysis_type_string=NULL;
     22        int      analysis_type;
     23        char*    sub_analysis_type_string=NULL;
     24        int      sub_analysis_type;
     25        int      numberofnodes;
    2526
    2627        /* output datasets: */
    … …  
    3940        FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
    4041        FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
     42        FetchData((void**)&numberofnodes,NULL,NULL,mxGetField(PARAMETERS,0,"numberofnodes"),"Integer",NULL);
    4143        FetchData((void**)&control_type,NULL,NULL,mxGetField(PARAMETERS,0,"control_type"),"String",NULL);
    4244        FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec");
    … …  
    5254
    5355        /*!Call core code: */
    54         Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
     56        Gradjx(&grad_g,numberofnodes,elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type);
    5557
    5658        /*write output : */
Note: See TracChangeset for help on using the changeset viewer.