Index: /issm/trunk/src/c/ControlConstrainx/ControlConstrainx.cpp =================================================================== --- /issm/trunk/src/c/ControlConstrainx/ControlConstrainx.cpp (revision 1045) +++ /issm/trunk/src/c/ControlConstrainx/ControlConstrainx.cpp (revision 1046) @@ -13,5 +13,5 @@ #include "../EnumDefinitions/EnumDefinitions.h" -void ControlConstrainx( double* p_g, int gsize, double mincontrolconstraint, double maxcontrolconstraint,char* control_type){ +void ControlConstrainx( double* p_g, int numdofnodes, double mincontrolconstraint, double maxcontrolconstraint,char* control_type){ int i; @@ -21,5 +21,5 @@ } else{ - for(i=0;iConfigure(elements,loads, nodes, materials); - /*Get size of matrix: */ - gsize=nodes->NumberOfDofs(); - /*Allocate grad_g: */ - grad_g=NewVec(gsize); + grad_g=NewVec(numberofnodes); /*Compute gradients: */ Index: /issm/trunk/src/c/Gradjx/Gradjx.h =================================================================== --- /issm/trunk/src/c/Gradjx/Gradjx.h (revision 1045) +++ /issm/trunk/src/c/Gradjx/Gradjx.h (revision 1046) @@ -9,5 +9,5 @@ /* local prototypes: */ -void Gradjx( Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, +void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials, double* u_g, double* lambda_g, ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type); Index: /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 1045) +++ /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 1046) @@ -135,10 +135,10 @@ parameters->AddObject(param); - param_g=(double*)xcalloc(model->numberofnodes*2,sizeof(double)); - for(i=0;inumberofnodes;i++)param_g[2*i+0]=control_parameter[i]; + param_g=(double*)xcalloc(model->numberofnodes,sizeof(double)); + for(i=0;inumberofnodes;i++)param_g[i]=control_parameter[i]; count++; param= new Param(count,"param_g",DOUBLEVEC); - param->SetDoubleVec(param_g,2*model->numberofnodes,2); + param->SetDoubleVec(param_g,model->numberofnodes,1); parameters->AddObject(param); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 1045) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 1046) @@ -4,5 +4,5 @@ #ifdef HAVE_CONFIG_H - #include "config.h" +#include "config.h" #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" @@ -33,9 +33,9 @@ Tria::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], - 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, - double tria_viscosity_overshoot,int tria_artdiff){ - + 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, + double tria_viscosity_overshoot,int tria_artdiff){ + int i; - + id=tria_id; mid=tria_mid; @@ -66,5 +66,5 @@ viscosity_overshoot=tria_viscosity_overshoot; artdiff=tria_artdiff; - + return; } @@ -156,8 +156,8 @@ /*get enum type of Tria: */ enum_type=TriaEnum(); - + /*marshall enum: */ memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); - + /*marshall Tria data: */ memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); @@ -187,37 +187,37 @@ memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot); memcpy(marshalled_dataset,&artdiff,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff); - + *pmarshalled_dataset=marshalled_dataset; return; } - + int Tria::MarshallSize(){ return sizeof(id) - +sizeof(mid) - +sizeof(mparid) - +sizeof(node_ids) - +sizeof(nodes) - +sizeof(node_offsets) - +sizeof(matice) - +sizeof(matice_offset) - +sizeof(matpar) - +sizeof(matpar_offset) - +sizeof(h) - +sizeof(s) - +sizeof(b) - +sizeof(k) - +sizeof(melting) - +sizeof(accumulation) - +sizeof(geothermalflux) - +sizeof(friction_type) - +sizeof(onbed) - +sizeof(p) - +sizeof(q) - +sizeof(shelf) - +sizeof(meanvel) - +sizeof(epsvel) - +sizeof(viscosity_overshoot) - +sizeof(artdiff) - +sizeof(int); //sizeof(int) for enum type + +sizeof(mid) + +sizeof(mparid) + +sizeof(node_ids) + +sizeof(nodes) + +sizeof(node_offsets) + +sizeof(matice) + +sizeof(matice_offset) + +sizeof(matpar) + +sizeof(matpar_offset) + +sizeof(h) + +sizeof(s) + +sizeof(b) + +sizeof(k) + +sizeof(melting) + +sizeof(accumulation) + +sizeof(geothermalflux) + +sizeof(friction_type) + +sizeof(onbed) + +sizeof(p) + +sizeof(q) + +sizeof(shelf) + +sizeof(meanvel) + +sizeof(epsvel) + +sizeof(viscosity_overshoot) + +sizeof(artdiff) + +sizeof(int); //sizeof(int) for enum type } @@ -291,5 +291,5 @@ int i; - + DataSet* loadsin=NULL; DataSet* nodesin=NULL; @@ -303,5 +303,5 @@ /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */ ResolvePointers((Object**)nodes,node_ids,node_offsets,3,nodesin); - + /*Same for materials: */ ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin); @@ -317,9 +317,9 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ if (analysis_type==ControlAnalysisEnum()){ - + CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==DiagnosticAnalysisEnum()){ - + if (sub_analysis_type==HorizAnalysisEnum()){ @@ -361,5 +361,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* gaussian points: */ int num_gauss,ig; @@ -375,5 +375,5 @@ double newviscosity; //viscosity double oldviscosity; //viscosity - + /* strain rate: */ double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ @@ -389,7 +389,7 @@ double Ke_gg[numdof][numdof]; //local element stiffness matrix double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - + double Jdet; - + /*input parameters for structural analysis (diagnostic): */ double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; @@ -414,10 +414,10 @@ for(i=0;iGetRhoIce()); printf(" gravity: %g \n",matpar->GetG()) - printf(" rho_water: %g \n",matpar->GetRhoWater()); + printf(" rho_water: %g \n",matpar->GetRhoWater()); printf(" Velocity: \n"); for (i=0;iGetViscosity2d(&viscosity, &epsilon[0]); matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - + /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ + onto this scalar matrix, so that we win some computational time: */ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=newviscosity*thickness*gauss_weight*Jdet; @@ -476,6 +476,6 @@ D[i][i]=D_scalar; } - - #ifdef _DEBUGELEMENTS_ + +#ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian loop %i\n",ig); @@ -492,5 +492,5 @@ printf(" gauss_weight: %g \n",gauss_weight); } - #endif +#endif /*Get B and Bprime matrices: */ @@ -500,12 +500,12 @@ /* Do the triple product tB*D*Bprime: */ TripleMultiply( &B[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof,0, - &Ke_gg_gaussian[0][0],0); + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof,0, + &Ke_gg_gaussian[0][0],0); /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ for( i=0; iRecover("dt",&dt); if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); @@ -645,5 +645,5 @@ v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); - + K[0][0]=pow(Jdettria,.5)/2.0*fabs(v_gauss[0]); K[1][1]=pow(Jdettria,.5)/2.0*fabs(v_gauss[1]); @@ -671,16 +671,16 @@ /* Do the triple product tL*D*L: */ TripleMultiply( &L[0],1,numdof,1, - &DL_scalar,1,1,0, - &L[0],1,numdof,0, - &Ke_gg_gaussian[0][0],0); - + &DL_scalar,1,1,0, + &L[0],1,numdof,0, + &Ke_gg_gaussian[0][0],0); + /*Get B and B prime matrix: */ GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); - + //Get vx, vy and their derivatives at gauss point GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); - + GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); @@ -702,12 +702,12 @@ TripleMultiply( &B[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &B[0][0],2,numdof,0, - &Ke_gg_thickness1[0][0],0); + &DL[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke_gg_thickness1[0][0],0); TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg_thickness2[0][0],0); + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_thickness2[0][0],0); /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ @@ -715,7 +715,7 @@ for( i=0; ielement_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); strcpy(friction->element_type,"2d"); - + friction->gravity=matpar->GetG(); friction->rho_ice=matpar->GetRhoIce(); @@ -974,14 +974,14 @@ DeleteFriction(&friction); - #ifdef _DEBUGELEMENTS_ +#ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); } - #endif +#endif /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - #ifdef _DEBUGELEMENTS_ +#ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian points: \n"); @@ -990,5 +990,5 @@ } } - #endif +#endif /* Start looping on the number of gaussian points: */ @@ -1025,5 +1025,5 @@ DL[i][i]=DL_scalar; } - + /* Do the triple producte tL*D*L: */ TripleMultiply( &L[0][0],2,numdof,1, @@ -1039,5 +1039,5 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - cleanup_and_return: +cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1062,5 +1062,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* gaussian points: */ int num_gauss,ig; @@ -1079,7 +1079,7 @@ double Ke_gg[numdof][numdof]; //local element stiffness matrix double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - + double Jdet; - + /* Get node coordinates and dof list: */ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); @@ -1100,5 +1100,5 @@ gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); - + /*Get L matrix: */ GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); @@ -1106,5 +1106,5 @@ /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - + DL_scalar=gauss_weight*Jdet; @@ -1121,6 +1121,6 @@ /*Add Ke_gg to global matrix Kgg: */ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - - cleanup_and_return: + +cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1132,21 +1132,21 @@ #define __FUNCT__ "Tria::CreatePVector" void Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ - + /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ if (analysis_type==ControlAnalysisEnum()){ - + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); - + } else if (analysis_type==DiagnosticAnalysisEnum()){ if (sub_analysis_type==HorizAnalysisEnum()){ - + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); - + } else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); } else if (analysis_type==SlopeComputeAnalysisEnum()){ - + CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); } @@ -1174,5 +1174,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* parameters: */ double plastic_stress; @@ -1215,5 +1215,5 @@ - #ifdef _DEBUGELEMENTS_ +#ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf("gravity %g\n",matpar->GetG()); @@ -1224,5 +1224,5 @@ printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]); } - #endif +#endif @@ -1230,5 +1230,5 @@ 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*/ - #ifdef _DEBUGELEMENTS_ +#ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian points: \n"); @@ -1237,5 +1237,5 @@ } } - #endif +#endif @@ -1251,7 +1251,7 @@ /*Compute thickness at gaussian point: */ GetParameterValue(&thickness, &h[0],gauss_l1l2l3); - + GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3); - + /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the * element itself: */ @@ -1262,6 +1262,6 @@ /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get nodal functions: */ + + /*Get nodal functions: */ GetNodalFunctions(l1l2l3, gauss_l1l2l3); @@ -1270,5 +1270,5 @@ - #ifdef _DEBUGELEMENTS_ +#ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian %i\n",ig); @@ -1280,5 +1280,5 @@ if(friction_type==1)printf(" plastic_stress(%g)\n",plastic_stress); } - #endif +#endif /*Build pe_g_gaussian vector: */ @@ -1303,5 +1303,5 @@ } //for (ig=0; igterms\n",ig); @@ -1315,10 +1315,10 @@ } } - #endif +#endif /*Add pe_g to global vector pg: */ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); - cleanup_and_return: +cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1342,5 +1342,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* gaussian points: */ int num_gauss,ig; @@ -1391,9 +1391,9 @@ GetParameterDerivativeValue(&slope[0], ¶m[0],&xyz_list[0][0], gauss_l1l2l3); - + /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get nodal functions: */ + + /*Get nodal functions: */ GetNodalFunctions(l1l2l3, gauss_l1l2l3); @@ -1414,5 +1414,5 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); - cleanup_and_return: +cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1445,5 +1445,5 @@ inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes); inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes); - + //Update material if necessary if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){ @@ -1452,5 +1452,5 @@ matice->SetB(B_average); } - + if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){ B_average=(B_list[0]+B_list[1]+B_list[2])/3.0; @@ -1459,5 +1459,5 @@ } - + void Tria::GetDofList(int* doflist,int* pnumberofdofspernode){ @@ -1465,5 +1465,5 @@ int doflist_per_node[MAXDOFSPERNODE]; int numberofdofspernode; - + for(i=0;i<3;i++){ nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); @@ -1491,8 +1491,8 @@ #define __FUNCT__ "Tria::GetParameterValue" void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){ - + /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian * point specifie by gauss_l1l2l3: */ - + /*nodal functions: */ double l1l2l3[3]; @@ -1513,5 +1513,5 @@ #define __FUNCT__ "Tria::GetParameterDerivativeValue" void Tria::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3){ - + const int NDOF2=2; const int numgrids=3; @@ -1523,5 +1523,5 @@ * p is a vector of size 2x1 already allocated. */ - + double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system. @@ -1548,5 +1548,5 @@ GetB(&B[0][0], xyz_list, gauss_l1l2l3); - #ifdef _DEBUG_ +#ifdef _DEBUG_ printf("B for grid1 : [ %lf %lf \n",B[0][0],B[0][1]); printf(" [ %lf %lf \n",B[1][0],B[1][1]); @@ -1560,14 +1560,14 @@ printf(" [ %lf %lf \n",B[1][4],B[1][5]); printf(" [ %lf %lf ]\n",B[2][4],B[2][5]); - + for (i=0;iparameters->FindParam((void*)&analysis_type,"analysis_type"); femmodel->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type"); + femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); femmodel->parameters->FindParam((void*)&solverstring,"solverstring"); femmodel->parameters->FindParam((void*)&control_type,"control_type"); @@ -69,5 +71,5 @@ VecToMPISerial(&lambda_g_double,lambda_g);VecFree(&lambda_g); - Gradjx( &grad_g, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, + Gradjx( &grad_g, numberofnodes,femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, u_g_double,lambda_g_double, inputs,analysis_type,sub_analysis_type,control_type); Index: /issm/trunk/src/c/parallel/ProcessResults.cpp =================================================================== --- /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 1045) +++ /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 1046) @@ -399,5 +399,5 @@ for(i=0;iAdd(control_type,param_g,2,numberofnodes); + inputs->Add(control_type,param_g,1,numberofnodes); inputs->Add("fit",fit[n]); @@ -123,9 +123,9 @@ _printf_("%s\n"," updating parameter using optimized search scalar..."); - for(i=0;iAdd(control_type,param_g,2,numberofnodes); + inputs->Add(control_type,param_g,1,numberofnodes); inputs->Add("fit",fit[n]); Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 1045) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 1046) @@ -64,15 +64,15 @@ /*First copy param_g so we don't modify it: */ - param_g_copy=(double*)xmalloc(gsize*sizeof(double)); - memcpy(param_g_copy,param_g,gsize*sizeof(double)); + param_g_copy=(double*)xmalloc(numberofnodes*sizeof(double)); + memcpy(param_g_copy,param_g,numberofnodes*sizeof(double)); /*First, update param_g using search_scalar: */ - for(i=0;iAdd(control_type,param_g_copy,2,numberofnodes); + inputs->Add(control_type,param_g_copy,1,numberofnodes); //Run diagnostic with updated parameters. @@ -85,5 +85,4 @@ u_g_double,u_g_obs, inputs,analysis_type,sub_analysis_type); - /*Free ressources:*/ xfree((void**)&optscal); Index: /issm/trunk/src/m/solutions/cielo/control_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/control_core.m (revision 1045) +++ /issm/trunk/src/m/solutions/cielo/control_core.m (revision 1046) @@ -18,5 +18,5 @@ %initialize control parameters, gradients and observations u_g_obs=m_dh.parameters.u_g_obs; -grad_g=zeros(m_dh.nodesets.gsize,1); +grad_g=zeros(m_dh.parameters.numberofnodes,1); param_g=models.dh.parameters.param_g; @@ -33,5 +33,5 @@ %update inputs with new fit inputs=add(inputs,'fit',m_dh.parameters.fit(n),'double'); - inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',2,m_dh.parameters.numberofnodes); + inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',1,m_dh.parameters.numberofnodes); %Update inputs in datasets Index: /issm/trunk/src/m/solutions/cielo/controlfinalsol.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/controlfinalsol.m (revision 1045) +++ /issm/trunk/src/m/solutions/cielo/controlfinalsol.m (revision 1046) @@ -6,5 +6,5 @@ %From parameters, build inputs for icediagnostic_core, using the final parameters -inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',2,m.parameters.numberofnodes); +inputs=add(inputs,m.parameters.control_type,param_g,'doublevec',1,m.parameters.numberofnodes); results.u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m (revision 1045) +++ /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m (revision 1046) @@ -10,5 +10,5 @@ %Plug parameter into inputs -inputs=add(inputs,m.parameters.control_type,parameter,'doublevec',2,m.parameters.numberofnodes); +inputs=add(inputs,m.parameters.control_type,parameter,'doublevec',1,m.parameters.numberofnodes); %Run diagnostic with updated parameters. Index: /issm/trunk/src/m/solutions/cielo/processresults.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/processresults.m (revision 1045) +++ /issm/trunk/src/m/solutions/cielo/processresults.m (revision 1046) @@ -124,5 +124,5 @@ newresults(i).step=results(i).step; newresults(i).time=results(i).time; - newresults(i).parameter=param_g(1:2:end); + newresults(i).parameter=param_g; else Index: /issm/trunk/src/mex/Gradj/Gradj.cpp =================================================================== --- /issm/trunk/src/mex/Gradj/Gradj.cpp (revision 1045) +++ /issm/trunk/src/mex/Gradj/Gradj.cpp (revision 1046) @@ -19,8 +19,9 @@ double* lambda_g=NULL; ParameterInputs* inputs=NULL; - char* analysis_type_string=NULL; - int analysis_type; - char* sub_analysis_type_string=NULL; - int sub_analysis_type; + char* analysis_type_string=NULL; + int analysis_type; + char* sub_analysis_type_string=NULL; + int sub_analysis_type; + int numberofnodes; /* output datasets: */ @@ -39,4 +40,5 @@ FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL); FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL); + FetchData((void**)&numberofnodes,NULL,NULL,mxGetField(PARAMETERS,0,"numberofnodes"),"Integer",NULL); FetchData((void**)&control_type,NULL,NULL,mxGetField(PARAMETERS,0,"control_type"),"String",NULL); FetchData((void**)&u_g,NULL,NULL,UG,"Vector","Vec"); @@ -52,5 +54,5 @@ /*!Call core code: */ - Gradjx(&grad_g, elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type); + Gradjx(&grad_g,numberofnodes,elements,nodes,loads,materials,u_g,lambda_g,inputs,analysis_type,sub_analysis_type,control_type); /*write output : */