Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 4046) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 4047) @@ -893,5 +893,5 @@ node=(Node*)(*object); - if (node->InAnalysis(analysis_type)){ + /*Check that this node corresponds to our analysis currently being carried out: */ /*Ok, this object is a node, ask it to plug values into partition: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 4046) +++ /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 4047) @@ -173,5 +173,5 @@ case PressureOldEnum : return "PressureOld"; case QmuPressureEnum : return "QmuPressure"; - case StokesPressureEnum : return "StokesPressure"; + case PressureStokesEnum : return "PressureStokes"; case ResetPenaltiesEnum : return "ResetPenalties"; case RheologyBEnum : return "RheologyB"; @@ -217,4 +217,5 @@ case CmMaxEnum : return "CmMax"; case CmMinEnum : return "CmMin"; + case GradientEnum : return "Gradient"; case ConnectivityEnum : return "Connectivity"; case ControlParameterEnum : return "ControlParameter"; Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4046) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 4047) @@ -252,4 +252,5 @@ CmMaxEnum, CmMinEnum, + GradientEnum, ConnectivityEnum, ControlParameterEnum, Index: /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 4046) +++ /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 4047) @@ -171,5 +171,5 @@ else if (strcmp(name,"PressureOld")==0) return PressureOldEnum; else if (strcmp(name,"QmuPressure")==0) return QmuPressureEnum; - else if (strcmp(name,"StokesPressure")==0) return StokesPressureEnum; + else if (strcmp(name,"PressureStokes")==0) return PressureStokesEnum; else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum; else if (strcmp(name,"RheologyB")==0) return RheologyBEnum; @@ -215,4 +215,5 @@ else if (strcmp(name,"CmMax")==0) return CmMaxEnum; else if (strcmp(name,"CmMin")==0) return CmMinEnum; + else if (strcmp(name,"Gradient")==0) return GradientEnum; else if (strcmp(name,"Connectivity")==0) return ConnectivityEnum; else if (strcmp(name,"ControlParameter")==0) return ControlParameterEnum; Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 4046) +++ /issm/trunk/src/c/Makefile.am (revision 4047) @@ -391,4 +391,6 @@ ./modules/DuplicateInputx/DuplicateInputx.h\ ./modules/DuplicateInputx/DuplicateInputx.cpp\ + ./modules/ScaleInputx/ScaleInputx.h\ + ./modules/ScaleInputx/ScaleInputx.cpp\ ./modules/ControlConstrainx/ControlConstrainx.h\ ./modules/ControlConstrainx/ControlConstrainx.cpp\ @@ -887,4 +889,6 @@ ./modules/DuplicateInputx/DuplicateInputx.h\ ./modules/DuplicateInputx/DuplicateInputx.cpp\ + ./modules/ScaleInputx/ScaleInputx.h\ + ./modules/ScaleInputx/ScaleInputx.cpp\ ./modules/ControlConstrainx/ControlConstrainx.h\ ./modules/ControlConstrainx/ControlConstrainx.cpp\ Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp =================================================================== --- /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp (revision 4046) +++ /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp (revision 4047) @@ -10,12 +10,9 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){ +void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){ - /*Intermediary*/ int i; - Element* element=NULL; + int dim; - /*output: */ - Vec grad_g=NULL; /*First, get elements and loads configured: */ @@ -24,18 +21,16 @@ parameters->Configure(elements,loads, nodes,vertices, materials,parameters); - /*Allocate grad_g: */ - grad_g=NewVec(numberofnodes); + /*retrieve some parameters: */ + parameters->FindParam(&dim,DimEnum); + /*Compute gradients: */ for (i=0;iSize();i++){ - element=(Element*)elements->GetObjectByOffset(i); - element->Gradj(grad_g,control_type); + Element* element=(Element*)elements->GetObjectByOffset(i); + element->Gradj(control_type); } - /*Assemble vector: */ - VecAssemblyBegin(grad_g); - VecAssemblyEnd(grad_g); + /*Extrude if needed: */ + if(dim==3) InputsExtrudex( elements,nodes,vertices,loads,materials,parameters,GradientEnum,0); - /*Assign output pointers: */ - *pgrad_g=grad_g; } Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.h =================================================================== --- /issm/trunk/src/c/modules/Gradjx/Gradjx.h (revision 4046) +++ /issm/trunk/src/c/modules/Gradjx/Gradjx.h (revision 4047) @@ -10,6 +10,5 @@ /* local prototypes: */ -void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, - int analysis_type,int sub_analysis_type,int control_type); +void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, int control_type); #endif /* _GRADJX_H */ Index: /issm/trunk/src/c/modules/ScaleInputx/ScaleInputx.cpp =================================================================== --- /issm/trunk/src/c/modules/ScaleInputx/ScaleInputx.cpp (revision 4047) +++ /issm/trunk/src/c/modules/ScaleInputx/ScaleInputx.cpp (revision 4047) @@ -0,0 +1,25 @@ +/*!\file ScaleInputx + * \brief: duplicte an input inside the elements, onto another, and wipe it off. + */ + +#include "./ScaleInputx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void ScaleInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int enum_type, double scale_factor){ + + /*intermediary:*/ + int i; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Go through elemnets, and ask to reinitialie the input: */ + for(i=0;iSize();i++){ + Element* element=(Element*)elements->GetObjectByOffset(i); + element->ScaleInput(enum_type,scale_factor); + } + +} Index: /issm/trunk/src/c/modules/ScaleInputx/ScaleInputx.h =================================================================== --- /issm/trunk/src/c/modules/ScaleInputx/ScaleInputx.h (revision 4047) +++ /issm/trunk/src/c/modules/ScaleInputx/ScaleInputx.h (revision 4047) @@ -0,0 +1,14 @@ +/*!\file: ScaleInputx.h + * \brief header file for field extrusion + */ + +#ifndef _SCALEINPUTX_H +#define _SCALEINPUTX_H + +#include "../../DataSet/DataSet.h" + +/* local prototypes: */ +void ScaleInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int enum_type, double scale_factor); + +#endif /* _SCALEINPUTX_H */ + Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 4046) +++ /issm/trunk/src/c/modules/modules.h (revision 4047) @@ -79,4 +79,5 @@ #include "./MaxAbsVzx/MaxAbsVzx.h" #include "./DuplicateInputx/DuplicateInputx.h" +#include "./ScaleInputx/ScaleInputx.h" #endif Index: /issm/trunk/src/c/objects/Elements/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4047) @@ -562,15 +562,15 @@ /*}}}*/ /*FUNCTION Beam::Gradj{{{1*/ -void Beam::Gradj(Vec, int control_type){ +void Beam::Gradj(int control_type){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Beam::GradjB{{{1*/ -void Beam::GradjB(Vec){ +void Beam::GradjB(void){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Beam::GradjDrag{{{1*/ -void Beam::GradjDrag(Vec){ +void Beam::GradjDrag(void){ ISSMERROR(" not supported yet!"); } @@ -990,2 +990,14 @@ } /*}}}*/ +/*FUNCTION Beam::ScaleInput(int enum_type,double scale_factor){{{1*/ +void Beam::ScaleInput(int enum_type,double scale_factor){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.h (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Beam.h (revision 4047) @@ -95,4 +95,5 @@ void MaxAbsVz(double* pmaxabsvz, bool process_units); void DuplicateInput(int original_enum,int new_enum); + void ScaleInput(int enum_type,double scale_factor); /*}}}*/ @@ -103,7 +104,7 @@ void GetThicknessList(double* thickness_list); void Du(Vec); - void Gradj(Vec, int control_type); - void GradjDrag(Vec); - void GradjB(Vec); + void Gradj(int control_type); + void GradjDrag(void); + void GradjB(void); double Misfit(void); double SurfaceArea(void); Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 4047) @@ -37,7 +37,7 @@ virtual void GetBedList(double* bed_list)=0; virtual void Du(Vec du_g)=0; - virtual void Gradj(Vec grad_g,int control_type)=0; - virtual void GradjDrag(Vec grad_g)=0; - virtual void GradjB(Vec grad_g)=0; + virtual void Gradj(int control_type)=0; + virtual void GradjDrag(void)=0; + virtual void GradjB(void)=0; virtual double Misfit(void)=0; virtual double CostFunction(void)=0; @@ -65,4 +65,5 @@ virtual void MaxAbsVz(double* pmaxabsvz, bool process_units)=0; virtual void DuplicateInput(int original_enum,int new_enum)=0; + virtual void ScaleInput(int enum_type,double scale_factor)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4047) @@ -4555,5 +4555,5 @@ /*}}}*/ /*FUNCTION Penta::Gradj {{{1*/ -void Penta::Gradj(Vec grad_g,int control_type){ +void Penta::Gradj(int control_type){ /*inputs: */ @@ -4567,8 +4567,8 @@ if (control_type==DragCoefficientEnum){ - GradjDrag( grad_g); + GradjDrag(); } else if (control_type=RheologyBEnum){ - GradjB( grad_g); + GradjB(); } else ISSMERROR("%s%i","control type not supported yet: ",control_type); @@ -4576,7 +4576,9 @@ /*}}}*/ /*FUNCTION Penta::GradjDrag {{{1*/ -void Penta::GradjDrag(Vec grad_g){ +void Penta::GradjDrag(void){ Tria* tria=NULL; + TriaVertexInput* triavertexinput=NULL; + double gradient[6]={0,0,0,0,0,0}; /*inputs: */ @@ -4586,5 +4588,5 @@ int analysis_type,sub_analysis_type; - /*retrive parameters: */ + /*retrieve parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); @@ -4608,7 +4610,5 @@ /*MacAyeal or Pattyn*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDrag( grad_g); - delete tria; - return; + tria->GradjDrag(); } else if (sub_analysis_type==StokesAnalysisEnum){ @@ -4616,15 +4616,23 @@ /*Stokes*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDragStokes( grad_g); - delete tria; - return; + tria->GradjDragStokes(); } else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); + + /*Now, the tria has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */ + triavertexinput=tria->inputs->GetInput(GradientEnum); + for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i]; + this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0])); + + delete tria; + } /*}}}*/ /*FUNCTION Penta::GradjB {{{1*/ -void Penta::GradjB(Vec grad_g){ +void Penta::GradjB(void){ Tria* tria=NULL; + TriaVertexInput* triavertexinput=NULL; + double gradient[6]={0,0,0,0,0,0}; /*inputs: */ @@ -4648,15 +4656,19 @@ * and compute gardj*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(grad_g); - delete tria; - return; + tria->GradjB(); } else{ /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(grad_g); - delete tria; - return; - } + tria->GradjB(); + } + + /*Now, the tria has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */ + triavertexinput=tria->inputs->GetInput(GradientEnum); + for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i]; + this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0])); + + delete tria; + } /*}}}*/ @@ -5328,2 +5340,14 @@ } /*}}}*/ +/*FUNCTION Penta::ScaleInput(int enum_type,double scale_factor){{{1*/ +void Penta::ScaleInput(int enum_type,double scale_factor){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 4047) @@ -86,7 +86,7 @@ bool GetOnBed(); void Du(Vec du_gg); - void Gradj(Vec grad_gg,int control_type); - void GradjDrag(Vec grad_gg); - void GradjB(Vec grad_gg); + void Gradj(int control_type); + void GradjDrag(void); + void GradjB(void); double Misfit(void); double SurfaceArea(void); @@ -162,4 +162,5 @@ void MaxAbsVz(double* pmaxabsvz, bool process_units); void DuplicateInput(int original_enum,int new_enum); + void ScaleInput(int enum_type,double scale_factor); Index: /issm/trunk/src/c/objects/Elements/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4047) @@ -383,15 +383,15 @@ /*}}}*/ /*FUNCTION Sing::Gradj {{{1*/ -void Sing::Gradj(Vec, int control_type){ +void Sing::Gradj(int control_type){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Sing::GradB {{{1*/ -void Sing::GradjB(Vec){ +void Sing::GradjB(void){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Sing::GradjDrag {{{1*/ -void Sing::GradjDrag(Vec){ +void Sing::GradjDrag(void){ ISSMERROR(" not supported yet!"); } @@ -696,2 +696,14 @@ } /*}}}*/ +/*FUNCTION Sing::ScaleInput(int enum_type,double scale_factor){{{1*/ +void Sing::ScaleInput(int enum_type,double scale_factor){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.h (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Sing.h (revision 4047) @@ -94,4 +94,5 @@ void MaxAbsVz(double* pmaxabsvz, bool process_units); void DuplicateInput(int original_enum,int new_enum); + void ScaleInput(int enum_type,double scale_factor); @@ -103,7 +104,7 @@ void GetThicknessList(double* thickness_list); void Du(Vec); - void Gradj(Vec, int control_type); - void GradjDrag(Vec); - void GradjB(Vec); + void Gradj(int control_type); + void GradjDrag(void); + void GradjB(void); double Misfit(void); double SurfaceArea(void); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4047) @@ -4125,5 +4125,5 @@ /*}}}*/ /*FUNCTION Tria::Gradj {{{1*/ -void Tria::Gradj(Vec grad_g,int control_type){ +void Tria::Gradj(int control_type){ /*inputs: */ @@ -4137,10 +4137,342 @@ if (control_type==DragCoefficientEnum){ - GradjDrag( grad_g); + GradjDrag(); } else if (control_type==RheologyBEnum){ - GradjB( grad_g); + GradjB(); } else ISSMERROR("%s%i","control type not supported yet: ",control_type); +} +/*}}}*/ +/*FUNCTION Tria::GradjDrag {{{1*/ +void Tria::GradjDrag(void){ + + + int i; + + /* node data: */ + const int numgrids=3; + const int NDOF2=2; + const int numdof=NDOF2*numgrids; + double xyz_list[numgrids][3]; + int doflist1[numgrids]; + double dh1dh3[NDOF2][numgrids]; + + /* grid data: */ + double adjx_list[numgrids]; + double adjy_list[numgrids]; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + /* parameters: */ + double dk[NDOF2]; + double vx,vy; + double lambda,mu; + double bed,thickness,Neff; + double alpha_complement; + int drag_type; + double drag; + Friction* friction=NULL; + + /*element vector at the gaussian points: */ + double grade_g[numgrids]={0.0}; + double grade_g_gaussian[numgrids]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3[3]; + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + + /*inputs: */ + bool shelf; + + /*parameters: */ + double cm_noisedmp; + double cm_mindmp_slope; + double cm_mindmp_value; + double cm_maxdmp_value; + double cm_maxdmp_slope; + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); + this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); + this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); + this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); + this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); + + + /*Get out if shelf*/ + if(shelf)return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList1(&doflist1[0]); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /* 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, 4); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); + else alpha_complement=0; + + /*Recover alpha_complement and k: */ + inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum); + + /*recover lambda and mu: */ + inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum); + inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum); + + /*recover vx and vy: */ + inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum); + inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Get nodal functions value at gaussian point:*/ + GetNodalFunctions(l1l2l3, gauss_l1l2l3); + + /*Get nodal functions derivatives*/ + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get k derivative: dk/dx */ + inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;icm_maxdmp_value){ + grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; + } + } + + /*Add gradje_g_gaussian vector to gradje_g: */ + for( i=0; iinputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + delete friction; + +} +/*}}}*/ +/*FUNCTION Tria::GradjDragStokes {{{1*/ +void Tria::GradjDragStokes(Vec grad_g){ + + int i; + + /* node data: */ + const int numgrids=3; + const int NDOF2=2; + double xyz_list[numgrids][3]; + int doflist1[numgrids]; + double dh1dh3[NDOF2][numgrids]; + + /* grid data: */ + double drag; + double alpha_complement; + Friction* friction=NULL; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + + /* parameters: */ + double vx,vy,vz; + double lambda,mu,xi; + double bed,thickness,Neff; + double surface_normal[3]; + double bed_normal[3]; + double dk[NDOF2]; + + /*element vector at the gaussian points: */ + double grade_g[numgrids]={0.0}; + double grade_g_gaussian[numgrids]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3[3]; + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + + /*inputs: */ + bool shelf; + int drag_type; + + /*parameters: */ + double cm_noisedmp; + double cm_mindmp_slope; + double cm_mindmp_value; + double cm_maxdmp_value; + double cm_maxdmp_slope; + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); + this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); + this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); + this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); + this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); + + /*Get out if shelf*/ + if(shelf)return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList1(&doflist1[0]); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /* 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, 4); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); + else alpha_complement=0; + inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum); + + /*recover lambda mu and xi: */ + inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum); + inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum); + inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum); + + /*recover vx vy and vz: */ + inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); + inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum); + + /*Get normal vector to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list); + + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result + bed_normal[1]=-surface_normal[1]; + bed_normal[2]=-surface_normal[2]; + + /* Get Jacobian determinant: */ + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Get nodal functions value at gaussian point:*/ + GetNodalFunctions(l1l2l3, gauss_l1l2l3); + + /*Get nodal functions derivatives*/ + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get k derivative: dk/dx */ + inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;icm_maxdmp_value){ + grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; + } + } + + /*Add gradje_g_gaussian vector to gradje_g: */ + for( i=0; iinputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + delete friction; + } /*}}}*/ @@ -4281,166 +4613,7 @@ } - /*Add grade_g to global vector grad_g: */ - VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES); - -cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); -} -/*}}}*/ -/*FUNCTION Tria::GradjDrag {{{1*/ -void Tria::GradjDrag(Vec grad_g){ - - - int i; - - /* node data: */ - const int numgrids=3; - const int NDOF2=2; - const int numdof=NDOF2*numgrids; - double xyz_list[numgrids][3]; - int doflist1[numgrids]; - double dh1dh3[NDOF2][numgrids]; - - /* grid data: */ - double adjx_list[numgrids]; - double adjy_list[numgrids]; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - /* parameters: */ - double dk[NDOF2]; - double vx,vy; - double lambda,mu; - double bed,thickness,Neff; - double alpha_complement; - int drag_type; - double drag; - Friction* friction=NULL; - - /*element vector at the gaussian points: */ - double grade_g[numgrids]={0.0}; - double grade_g_gaussian[numgrids]; - - /* Jacobian: */ - double Jdet; - - /*nodal functions: */ - double l1l2l3[3]; - - /* strain rate: */ - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - - /*inputs: */ - bool shelf; - - /*parameters: */ - double cm_noisedmp; - double cm_mindmp_slope; - double cm_mindmp_value; - double cm_maxdmp_value; - double cm_maxdmp_slope; - - int analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); - this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); - this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); - this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); - - - /*Get out if shelf*/ - if(shelf)return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList1(&doflist1[0]); - - /*Build frictoin element, needed later: */ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - friction=new Friction("2d",inputs,matpar,analysis_type); - - /* 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, 4); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); - else alpha_complement=0; - - /*Recover alpha_complement and k: */ - inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum); - - /*recover lambda and mu: */ - inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum); - inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum); - - /*recover vx and vy: */ - inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum); - inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - - /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); - - /*Get k derivative: dk/dx */ - inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;icm_maxdmp_value){ - grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; - } - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for( i=0; iinputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); + cleanup_and_return: xfree((void**)&first_gauss_area_coord); @@ -4448,176 +4621,4 @@ xfree((void**)&third_gauss_area_coord); xfree((void**)&gauss_weights); - delete friction; - -} -/*}}}*/ -/*FUNCTION Tria::GradjDragStokes {{{1*/ -void Tria::GradjDragStokes(Vec grad_g){ - - int i; - - /* node data: */ - const int numgrids=3; - const int NDOF2=2; - double xyz_list[numgrids][3]; - int doflist1[numgrids]; - double dh1dh3[NDOF2][numgrids]; - - /* grid data: */ - double drag; - double alpha_complement; - Friction* friction=NULL; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - - /* parameters: */ - double vx,vy,vz; - double lambda,mu,xi; - double bed,thickness,Neff; - double surface_normal[3]; - double bed_normal[3]; - double dk[NDOF2]; - - /*element vector at the gaussian points: */ - double grade_g[numgrids]={0.0}; - double grade_g_gaussian[numgrids]; - - /* Jacobian: */ - double Jdet; - - /*nodal functions: */ - double l1l2l3[3]; - - /* strain rate: */ - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - - /*inputs: */ - bool shelf; - int drag_type; - - /*parameters: */ - double cm_noisedmp; - double cm_mindmp_slope; - double cm_mindmp_value; - double cm_maxdmp_value; - double cm_maxdmp_slope; - - int analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - inputs->GetParameterValue(&drag_type,DragTypeEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); - this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); - this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); - this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); - - /*Get out if shelf*/ - if(shelf)return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList1(&doflist1[0]); - - /*Build frictoin element, needed later: */ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - friction=new Friction("2d",inputs,matpar,analysis_type); - - /* 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, 4); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); - else alpha_complement=0; - inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum); - - /*recover lambda mu and xi: */ - inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum); - inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum); - inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum); - - /*recover vx vy and vz: */ - inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); - inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum); - - /*Get normal vector to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list); - - bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result - bed_normal[1]=-surface_normal[1]; - bed_normal[2]=-surface_normal[2]; - - /* Get Jacobian determinant: */ - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - - /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); - - /*Get k derivative: dk/dx */ - inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;icm_maxdmp_value){ - grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; - } - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for( i=0; iinputs->GetInput(enum_type); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 4046) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 4047) @@ -88,9 +88,9 @@ void GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3); void Du(Vec du_g); - void Gradj(Vec grad_g,int control_type); - void GradjDrag(Vec grad_g); - void GradjDragStokes(Vec grad_g); + void Gradj(int control_type); + void GradjDrag(void); + void GradjDragStokes(void); + void GradjB(void); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); - void GradjB(Vec grad_g); double Misfit(void); double SurfaceArea(void); @@ -140,4 +140,5 @@ void MaxAbsVz(double* pmaxabsvz, bool process_units); void DuplicateInput(int original_enum,int new_enum); + void ScaleInput(int enum_type,double scale_factor); Index: /issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp (revision 4046) +++ /issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp (revision 4047) @@ -250,2 +250,11 @@ } /*}}}*/ +/*FUNCTION BeamVertexInput::Scale(double scale_factor){{{1*/ +void BeamVertexInput::Scale(double scale_factor){ + + int i; + const int numgrids=2; + + for(i=0;iparameters->FindParam(&solverstring,SolverStringEnum); + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&isstokes,IsStokesEnum); + + /*set analysis type to compute velocity: */ + if(isstokes)femmodel->SetAnalysis(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum); + else femmodel->SetAnalysis(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum); + + _printf_("%s\n"," recover solution for this stiffness and right hand side:"); + solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0,NULL, femmodel,DiagnosticAnalysisEnum,sub_analysis_type); + + _printf_("%s\n"," buid Du, difference between observed velocity and model velocity:"); + Dux(du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + _printf_("%s\n"," reduce adjoint load from g-set to f-set:"); + Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys, femmodel->nodesets,true); //true means that ys0 flag is activated: all spcs show 0 displacement + + /*free some ressources: */ + VecFree(&du_g);MatFree(&K_fs0); + + _printf_("%s\n"," solve for adjoint vector:"); + Solverx(&adjoint_f, K_ff0, du_f, NULL, solverstring); + + /*free some ressources: */ + VecFree(&du_f); MatFree(&K_ff0); + + _printf_("%s\n"," merge back to g set:"); + Mergesolutionfromftogx(&adjoint_g, adjoint_f,femmodel->Gmn,femmodel->ys,femmodel->nodesets,true);//true means that ys0 flag is activated: all spc are 0 + + /*free some ressources: */ + VecFree(&adjoint_f); + + /*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */ + if(isstokes)femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum); + else femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum); + + UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g,AdjointAnalysisEnum, NoneAnalysisEnum); + + /*Free ressources:*/ + char* solverstring=NULL; + VecFree(&ug); + MatFree(&K_ff0); + MatFree(&Mat K_fs0); + VecFree(&du_g); + VecFree(&du_f); + VecFree(&adjoint_f); + VecFree(&adjoint_g); + +} Index: /issm/trunk/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/control_core.cpp (revision 4046) +++ /issm/trunk/src/c/solutions/control_core.cpp (revision 4047) @@ -3,23 +3,15 @@ */ +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" #include "./solutions.h" #include "../modules/modules.h" -#include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/include.h" #include "../solvers/solvers.h" -Results* control_core(Model* model){ - - extern int my_rank; - - /*fem model: */ - FemModel* fem_model=NULL; - - /*output: */ - Results* results=NULL; - - /*Intermediary: */ - Results* diagnostic_results=NULL; - Results* gradjcompute_results=NULL; - Results* steadystate_results=NULL; +void control_core(FemModel* femmodel){ + Vec u_g=NULL; Vec t_g=NULL; @@ -29,5 +21,4 @@ double* fit=NULL; double* optscal=NULL; - int gsize; double* maxiter=NULL; double* cm_jump=NULL; @@ -56,29 +47,24 @@ int numberofnodes; - //initialize results - results=new Results(); - /*Process models*/ - ControlInitialization(model); - fem_model=model->GetActiveFormulation(); + ControlInitialization(femmodel); /*Recover parameters used throughout the solution:*/ - model->FindParam(&nsteps,NStepsEnum); - model->FindParam(&control_type,ControlTypeEnum); - model->FindParam(&fit,NULL,NULL,FitEnum); - model->FindParam(&optscal,NULL,NULL,OptScalEnum); - model->FindParam(&maxiter,NULL,NULL,MaxIterEnum); - model->FindParam(&cm_jump,NULL,NULL,CmJumpEnum); - model->FindParam(&eps_cm,EpsCmEnum); - model->FindParam(&tolx,TolXEnum); - model->FindParam(&cm_min,CmMinEnum); - model->FindParam(&cm_max,CmMaxEnum); - model->FindParam(&cm_gradient,CmGradientEnum); - model->FindParam(¶m_g,NULL,NULL,ControlParameterEnum); - model->FindParam(&analysis_type,AnalysisTypeEnum); - model->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&control_steady,ControlSteadyEnum); - gsize=fem_model->nodes->NumberOfDofs(); + femmodel->parameters->FindParam(&nsteps,NStepsEnum); + femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&fit,NULL,NULL,FitEnum); + femmodel->parameters->FindParam(&optscal,NULL,NULL,OptScalEnum); + femmodel->parameters->FindParam(&maxiter,NULL,NULL,MaxIterEnum); + femmodel->parameters->FindParam(&cm_jump,NULL,NULL,CmJumpEnum); + femmodel->parameters->FindParam(&eps_cm,EpsCmEnum); + femmodel->parameters->FindParam(&tolx,TolXEnum); + femmodel->parameters->FindParam(&cm_min,CmMinEnum); + femmodel->parameters->FindParam(&cm_max,CmMaxEnum); + femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum); + femmodel->parameters->FindParam(¶m_g,NULL,NULL,ControlParameterEnum); + femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); + femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); /*Initialize misfit: */ @@ -89,41 +75,19 @@ _printf_("\n%s%i%s%i\n"," control method step ",n+1,"/",nsteps); - model->UpdateInputsFromVector(param_g,control_type,VertexEnum); - model->UpdateInputsFromConstant(fit[n],FitEnum); + inputs->Add(control_type,param_g,1,numberofnodes); + femmodel->UpdateInputsFromConstant(fit[n],FitEnum); /*In case we are running a steady state control method, compute new temperature field using new parameter * distribution: */ - if (control_steady){ - steadystate_results= steadystate_core(model); - VecFree(&t_g); steadystate_results->FindResult(&t_g,"t_g"); - delete steadystate_results; - model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum); - } + if (control_steady) steadystate_core(model); - _printf_("%s\n"," computing gradJ..."); - gradjcompute_results= gradjcompute_core(model); - gradjcompute_results->FindResult(&grad_g,"grad_g"); - delete gradjcompute_results; + if(verbose)_printf_("%s\n"," computing gradJ..."); + gradient_core(femmodel); /*Return gradient if asked: */ if (cm_gradient){ - - /*Plug results into output dataset: */ - results->AddObject(new Result(results->Size()+1,0,1,"grad_g",grad_g)); - results->AddObject(new Result(results->Size()+1,0,1,"analysis_type",EnumAsString(analysis_type))); - - /*Free ressources: */ - xfree((void**)&control_type); - xfree((void**)&fit); - xfree((void**)&optscal); - xfree((void**)&maxiter); - xfree((void**)&cm_jump); - xfree((void**)&grad_g_double); - xfree((void**)¶m_g); - VecFree(&u_g); - VecFree(&t_g); - VecFree(&m_g); - xfree((void**)&J); - return results; + /*Transfer gradient from input to results: */ + InputToResultx(femodel->elements,femodel->nodes,femodel->vertices,femodel->loads,femodel->materials,femodel->parameters,GradientEnum); + return; } Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4046) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4047) @@ -68,5 +68,5 @@ //"recondition" pressure computed previously: DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); - ScaleInputx(femmmodel,PressureStokesEnum,1.0/stokesreconditioning); + ScaleInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning); if(verbose)_printf_("%s\n"," update boundary conditions for stokes using velocities previously computed..."); @@ -74,5 +74,5 @@ if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ..."); - solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,fem_ds,DiagnosticAnalysisEnum,StokesAnalysisEnum); + solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum); } } Index: /issm/trunk/src/c/solutions/gradient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/gradient_core.cpp (revision 4047) +++ /issm/trunk/src/c/solutions/gradient_core.cpp (revision 4047) @@ -0,0 +1,36 @@ +/*!\file: gradient_core.cpp + * \brief compute inverse method gradient direction. + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +#include "../solvers/solvers.h" + + +void gradient_core(FemModel* femmodel){ + + + /*parameters: */ + int control_steady; + int control_type; + int verbose; + + /*retrieve parameters:*/ + femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); + femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&verbose,VerboseEnum); + + if(verbose)_printf_("%s\n"," compute adjoint state:"); + adjoint_core(femmodel); + + if(verbose)_printf_("%s\n"," compute gradient:"); + Gradjx( femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type); + + if(control_steady)diagnostic_core(model); + +} Index: sm/trunk/src/c/solutions/gradjcompute_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/gradjcompute_core.cpp (revision 4046) +++ (revision ) @@ -1,135 +1,0 @@ -/*!\file: gradjcompute_core.cpp - * \brief compute inverse method gradient direction. - */ - -#include "../modules/modules.h" -#include "./solutions.h" -#include "../EnumDefinitions/EnumDefinitions.h" - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "../solvers/solvers.h" - -Results* gradjcompute_core(Model* model){ - - - /*intermediary: */ - FemModel* femmodel=NULL; - Results* diagnostic_results=NULL; - Result* result=NULL; - int analysis_type; - int sub_analysis_type; - int numberofnodes; - int control_steady; - int numberofdofspernode; - char* solverstring=NULL; - int control_type; - - Vec u_g=NULL; - - Vec du_g=NULL; - Vec du_f=NULL; - Vec grad_g=NULL; - - Vec lambda_f=NULL; - Vec lambda_g=NULL; - double* lambdax=NULL; - double* lambday=NULL; - double* vx=NULL; - double* vy=NULL; - double* vz=NULL; - - Mat K_ff0=NULL; - Mat K_fs0=NULL; - - /*output: */ - Results* results=NULL; - - /*flags: */ - int verbose=0; - int dim=-1; - int extrude_param=0; - - //initialize results - results=new Results(); - - /*some parameters:*/ - femmodel=model->GetActiveFormulation(); - femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum); - femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); - femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); - femmodel->parameters->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum); - femmodel->parameters->FindParam(&solverstring,SolverStringEnum); - femmodel->parameters->FindParam(&control_type,ControlTypeEnum); - femmodel->parameters->FindParam(&extrude_param,ExtrudeParamEnum); - femmodel->parameters->FindParam(&verbose,VerboseEnum); - femmodel->parameters->FindParam(&dim,DimEnum); - - _printf_("%s\n"," recover solution for this stiffness and right hand side:"); - diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,NULL, femmodel,DiagnosticAnalysisEnum,sub_analysis_type); - - _printf_("%s\n"," buid Du, difference between observed velocity and model velocity:"); - Dux( &du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,analysis_type,sub_analysis_type); - - _printf_("%s\n"," reduce adjoint load from g-set to f-set:"); - Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys, femmodel->nodesets,true); //true means that ys0 flag is activated: all spc are 0 - VecFree(&du_g);MatFree(&K_fs0); - - _printf_("%s\n"," solve for adjoint vector:"); - Solverx(&lambda_f, K_ff0, du_f, NULL, solverstring); - VecFree(&du_f); MatFree(&K_ff0); - - _printf_("%s\n"," merge back to g set:"); - Mergesolutionfromftogx(&lambda_g, lambda_f,femmodel->Gmn,femmodel->ys,femmodel->nodesets,true);//true means that ys0 flag is activated: all spc are 0 - VecFree(&lambda_f); - - /*add to inputs: */ - SplitSolutionVectorx(lambda_g,numberofnodes,numberofdofspernode,&lambdax,&lambday); - model->UpdateInputsFromVector(lambdax,AdjointxEnum,VertexEnum); - model->UpdateInputsFromVector(lambday,AdjointyEnum,VertexEnum); - VecFree(&lambda_g); - - _printf_("%s\n"," compute gradJ:"); - Gradjx( &grad_g, numberofnodes,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, - analysis_type,sub_analysis_type,control_type); - - if (dim==3 && extrude_param){ - - _printf_("%s\n"," extruding gradient..."); - FieldExtrudex( grad_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"gradj",0); - } - - if(control_steady){ - diagnostic_results= diagnostic_core(model); - - //extract u_g and add it to input (3d velocity needed by thermal_core) - diagnostic_results->FindResult(&u_g,"u_g"); - - SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz); - model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); - model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); - - delete diagnostic_results; - } - - /*Plug results into output dataset: */ - InputToResultx(&result,femmodel->elements,femmodel->nodes,femmodel->vertices, femmodel->loads, femmodel->materials,femmodel->parameters,GradientAnalysisEnum,results->Size()+1,0,1); results->AddObject(result); - results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(GradientAnalysisEnum))); - - /*Free ressources:*/ - VecFree(&u_g); - VecFree(&grad_g); - xfree((void**)&solverstring); - xfree((void**)&lambdax); - xfree((void**)&lambday); - xfree((void**)&vx); - xfree((void**)&vy); - xfree((void**)&vz); - -} Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 4046) +++ /issm/trunk/src/c/solutions/solutions.h (revision 4047) @@ -13,6 +13,4 @@ /*cores: */ -Results* gradjcompute_core(FemModel* model); -void diagnostic_core(FemModel* model); Results* prognostic_core(FemModel* model); Results* prognostic2_core(FemModel* model); @@ -26,5 +24,8 @@ Results* transient_core_2d(FemModel* model); Results* transient_core_3d(FemModel* model); -Results* thermal_core(FemModel* model); +void adjoint_core(FemModel* model); +void gradient_core(FemModel* model); +void diagnostic_core(FemModel* model); +void thermal_core(FemModel* model); void surfaceslope_core(FemModel* femmodel); void bedslope_core(FemModel* femmodel); Index: /issm/trunk/src/c/solvers/solvers.h =================================================================== --- /issm/trunk/src/c/solvers/solvers.h (revision 4046) +++ /issm/trunk/src/c/solvers/solvers.h (revision 4047) @@ -15,4 +15,5 @@ void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type); void solver_linear(Vec* pug, FemModel* femmodel,int analysis_type,int sub_analysis_type); +void solver_adjoint(Vec* pug, FemModel* femmodel,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/m/enum/ConnectivityEnum.m =================================================================== --- /issm/trunk/src/m/enum/ConnectivityEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/ConnectivityEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=ConnectivityEnum() -macro=201; +macro=202; Index: /issm/trunk/src/m/enum/ControlParameterEnum.m =================================================================== --- /issm/trunk/src/m/enum/ControlParameterEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/ControlParameterEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=ControlParameterEnum() -macro=202; +macro=203; Index: /issm/trunk/src/m/enum/ControlSteadyEnum.m =================================================================== --- /issm/trunk/src/m/enum/ControlSteadyEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/ControlSteadyEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=ControlSteadyEnum() -macro=203; +macro=204; Index: /issm/trunk/src/m/enum/DakotaParameterEnum.m =================================================================== --- /issm/trunk/src/m/enum/DakotaParameterEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/DakotaParameterEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=DakotaParameterEnum() -macro=204; +macro=205; Index: /issm/trunk/src/m/enum/DimEnum.m =================================================================== --- /issm/trunk/src/m/enum/DimEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/DimEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=DimEnum() -macro=205; +macro=206; Index: /issm/trunk/src/m/enum/EnumAsString.m =================================================================== --- /issm/trunk/src/m/enum/EnumAsString.m (revision 4046) +++ /issm/trunk/src/m/enum/EnumAsString.m (revision 4047) @@ -168,5 +168,5 @@ case PressureOldEnum(), string='PressureOld'; return case QmuPressureEnum(), string='QmuPressure'; return - case StokesPressureEnum(), string='StokesPressure'; return + case PressureStokesEnum(), string='PressureStokes'; return case ResetPenaltiesEnum(), string='ResetPenalties'; return case RheologyBEnum(), string='RheologyB'; return @@ -212,4 +212,5 @@ case CmMaxEnum(), string='CmMax'; return case CmMinEnum(), string='CmMin'; return + case GradientEnum(), string='Gradient'; return case ConnectivityEnum(), string='Connectivity'; return case ControlParameterEnum(), string='ControlParameter'; return Index: /issm/trunk/src/m/enum/EpsAbsEnum.m =================================================================== --- /issm/trunk/src/m/enum/EpsAbsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/EpsAbsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=EpsAbsEnum() -macro=206; +macro=207; Index: /issm/trunk/src/m/enum/EpsCmEnum.m =================================================================== --- /issm/trunk/src/m/enum/EpsCmEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/EpsCmEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=EpsCmEnum() -macro=207; +macro=208; Index: /issm/trunk/src/m/enum/EpsRelEnum.m =================================================================== --- /issm/trunk/src/m/enum/EpsRelEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/EpsRelEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=EpsRelEnum() -macro=208; +macro=209; Index: /issm/trunk/src/m/enum/EpsResEnum.m =================================================================== --- /issm/trunk/src/m/enum/EpsResEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/EpsResEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=EpsResEnum() -macro=209; +macro=210; Index: /issm/trunk/src/m/enum/ExtrudeParamEnum.m =================================================================== --- /issm/trunk/src/m/enum/ExtrudeParamEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/ExtrudeParamEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=ExtrudeParamEnum() -macro=210; +macro=211; Index: /issm/trunk/src/m/enum/GradientEnum.m =================================================================== --- /issm/trunk/src/m/enum/GradientEnum.m (revision 4047) +++ /issm/trunk/src/m/enum/GradientEnum.m (revision 4047) @@ -0,0 +1,11 @@ +function macro=GradientEnum() +%GRADIENTENUM - Enum of Gradient +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/SynchronizeMatlabEnum +% Please read src/c/README for more information +% +% Usage: +% macro=GradientEnum() + +macro=201; Index: /issm/trunk/src/m/enum/HeatCapacityEnum.m =================================================================== --- /issm/trunk/src/m/enum/HeatCapacityEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/HeatCapacityEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=HeatCapacityEnum() -macro=211; +macro=212; Index: /issm/trunk/src/m/enum/IsHutterEnum.m =================================================================== --- /issm/trunk/src/m/enum/IsHutterEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/IsHutterEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=IsHutterEnum() -macro=212; +macro=213; Index: /issm/trunk/src/m/enum/IsMacAyealPattynEnum.m =================================================================== --- /issm/trunk/src/m/enum/IsMacAyealPattynEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/IsMacAyealPattynEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=IsMacAyealPattynEnum() -macro=213; +macro=214; Index: /issm/trunk/src/m/enum/IsStokesEnum.m =================================================================== --- /issm/trunk/src/m/enum/IsStokesEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/IsStokesEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=IsStokesEnum() -macro=214; +macro=215; Index: /issm/trunk/src/m/enum/LatentHeatEnum.m =================================================================== --- /issm/trunk/src/m/enum/LatentHeatEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/LatentHeatEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=LatentHeatEnum() -macro=215; +macro=216; Index: /issm/trunk/src/m/enum/LowmemEnum.m =================================================================== --- /issm/trunk/src/m/enum/LowmemEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/LowmemEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=LowmemEnum() -macro=216; +macro=217; Index: /issm/trunk/src/m/enum/MaxIterEnum.m =================================================================== --- /issm/trunk/src/m/enum/MaxIterEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/MaxIterEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=MaxIterEnum() -macro=217; +macro=218; Index: /issm/trunk/src/m/enum/MaxNonlinearIterationsEnum.m =================================================================== --- /issm/trunk/src/m/enum/MaxNonlinearIterationsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/MaxNonlinearIterationsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=MaxNonlinearIterationsEnum() -macro=218; +macro=219; Index: /issm/trunk/src/m/enum/MeltingPointEnum.m =================================================================== --- /issm/trunk/src/m/enum/MeltingPointEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/MeltingPointEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=MeltingPointEnum() -macro=219; +macro=220; Index: /issm/trunk/src/m/enum/MinMechanicalConstraintsEnum.m =================================================================== --- /issm/trunk/src/m/enum/MinMechanicalConstraintsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/MinMechanicalConstraintsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=MinMechanicalConstraintsEnum() -macro=220; +macro=221; Index: /issm/trunk/src/m/enum/MinThermalConstraintsEnum.m =================================================================== --- /issm/trunk/src/m/enum/MinThermalConstraintsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/MinThermalConstraintsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=MinThermalConstraintsEnum() -macro=221; +macro=222; Index: /issm/trunk/src/m/enum/NStepsEnum.m =================================================================== --- /issm/trunk/src/m/enum/NStepsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NStepsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NStepsEnum() -macro=222; +macro=223; Index: /issm/trunk/src/m/enum/NdtEnum.m =================================================================== --- /issm/trunk/src/m/enum/NdtEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NdtEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NdtEnum() -macro=223; +macro=224; Index: /issm/trunk/src/m/enum/NumOutputEnum.m =================================================================== --- /issm/trunk/src/m/enum/NumOutputEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NumOutputEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NumOutputEnum() -macro=224; +macro=225; Index: /issm/trunk/src/m/enum/NumRiftsEnum.m =================================================================== --- /issm/trunk/src/m/enum/NumRiftsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NumRiftsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NumRiftsEnum() -macro=225; +macro=226; Index: /issm/trunk/src/m/enum/NumberOfDofsPerNodeEnum.m =================================================================== --- /issm/trunk/src/m/enum/NumberOfDofsPerNodeEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NumberOfDofsPerNodeEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NumberOfDofsPerNodeEnum() -macro=226; +macro=227; Index: /issm/trunk/src/m/enum/NumberOfElementsEnum.m =================================================================== --- /issm/trunk/src/m/enum/NumberOfElementsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NumberOfElementsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NumberOfElementsEnum() -macro=227; +macro=228; Index: /issm/trunk/src/m/enum/NumberOfNodesEnum.m =================================================================== --- /issm/trunk/src/m/enum/NumberOfNodesEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NumberOfNodesEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NumberOfNodesEnum() -macro=228; +macro=229; Index: /issm/trunk/src/m/enum/NumberOfVerticesEnum.m =================================================================== --- /issm/trunk/src/m/enum/NumberOfVerticesEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/NumberOfVerticesEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=NumberOfVerticesEnum() -macro=229; +macro=230; Index: /issm/trunk/src/m/enum/OptScalEnum.m =================================================================== --- /issm/trunk/src/m/enum/OptScalEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/OptScalEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=OptScalEnum() -macro=230; +macro=231; Index: /issm/trunk/src/m/enum/OutputFileNameEnum.m =================================================================== --- /issm/trunk/src/m/enum/OutputFileNameEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/OutputFileNameEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=OutputFileNameEnum() -macro=231; +macro=232; Index: /issm/trunk/src/m/enum/ParameterOutputEnum.m =================================================================== --- /issm/trunk/src/m/enum/ParameterOutputEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/ParameterOutputEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=ParameterOutputEnum() -macro=232; +macro=233; Index: /issm/trunk/src/m/enum/PenaltyMeltingEnum.m =================================================================== --- /issm/trunk/src/m/enum/PenaltyMeltingEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/PenaltyMeltingEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=PenaltyMeltingEnum() -macro=233; +macro=234; Index: /issm/trunk/src/m/enum/PressureStokesEnum.m =================================================================== --- /issm/trunk/src/m/enum/PressureStokesEnum.m (revision 4047) +++ /issm/trunk/src/m/enum/PressureStokesEnum.m (revision 4047) @@ -0,0 +1,11 @@ +function macro=PressureStokesEnum() +%PRESSURESTOKESENUM - Enum of PressureStokes +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/SynchronizeMatlabEnum +% Please read src/c/README for more information +% +% Usage: +% macro=PressureStokesEnum() + +macro=157; Index: /issm/trunk/src/m/enum/QmuAnalysisEnum.m =================================================================== --- /issm/trunk/src/m/enum/QmuAnalysisEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/QmuAnalysisEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=QmuAnalysisEnum() -macro=234; +macro=235; Index: /issm/trunk/src/m/enum/QmuErrNameEnum.m =================================================================== --- /issm/trunk/src/m/enum/QmuErrNameEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/QmuErrNameEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=QmuErrNameEnum() -macro=235; +macro=236; Index: /issm/trunk/src/m/enum/QmuInNameEnum.m =================================================================== --- /issm/trunk/src/m/enum/QmuInNameEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/QmuInNameEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=QmuInNameEnum() -macro=236; +macro=237; Index: /issm/trunk/src/m/enum/QmuMassFluxSegmentsEnum.m =================================================================== --- /issm/trunk/src/m/enum/QmuMassFluxSegmentsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/QmuMassFluxSegmentsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=QmuMassFluxSegmentsEnum() -macro=237; +macro=238; Index: /issm/trunk/src/m/enum/QmuNPartEnum.m =================================================================== --- /issm/trunk/src/m/enum/QmuNPartEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/QmuNPartEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=QmuNPartEnum() -macro=238; +macro=239; Index: /issm/trunk/src/m/enum/QmuOutNameEnum.m =================================================================== --- /issm/trunk/src/m/enum/QmuOutNameEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/QmuOutNameEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=QmuOutNameEnum() -macro=239; +macro=240; Index: /issm/trunk/src/m/enum/QmuPartEnum.m =================================================================== --- /issm/trunk/src/m/enum/QmuPartEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/QmuPartEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=QmuPartEnum() -macro=240; +macro=241; Index: /issm/trunk/src/m/enum/ResponseDescriptorsEnum.m =================================================================== --- /issm/trunk/src/m/enum/ResponseDescriptorsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/ResponseDescriptorsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=ResponseDescriptorsEnum() -macro=241; +macro=242; Index: /issm/trunk/src/m/enum/SolverStringEnum.m =================================================================== --- /issm/trunk/src/m/enum/SolverStringEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/SolverStringEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=SolverStringEnum() -macro=242; +macro=243; Index: /issm/trunk/src/m/enum/SparsityEnum.m =================================================================== --- /issm/trunk/src/m/enum/SparsityEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/SparsityEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=SparsityEnum() -macro=243; +macro=244; Index: sm/trunk/src/m/enum/StokesPressureEnum.m =================================================================== --- /issm/trunk/src/m/enum/StokesPressureEnum.m (revision 4046) +++ (revision ) @@ -1,11 +1,0 @@ -function macro=StokesPressureEnum() -%STOKESPRESSUREENUM - Enum of StokesPressure -% -% WARNING: DO NOT MODIFY THIS FILE -% this file has been automatically generated by src/c/SynchronizeMatlabEnum -% Please read src/c/README for more information -% -% Usage: -% macro=StokesPressureEnum() - -macro=157; Index: /issm/trunk/src/m/enum/StringAsEnum.m =================================================================== --- /issm/trunk/src/m/enum/StringAsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/StringAsEnum.m (revision 4047) @@ -166,5 +166,5 @@ elseif (strcmpi(name,'PressureOld')), enum=PressureOldEnum(); return elseif (strcmpi(name,'QmuPressure')), enum=QmuPressureEnum(); return -elseif (strcmpi(name,'StokesPressure')), enum=StokesPressureEnum(); return +elseif (strcmpi(name,'PressureStokes')), enum=PressureStokesEnum(); return elseif (strcmpi(name,'ResetPenalties')), enum=ResetPenaltiesEnum(); return elseif (strcmpi(name,'RheologyB')), enum=RheologyBEnum(); return @@ -210,4 +210,5 @@ elseif (strcmpi(name,'CmMax')), enum=CmMaxEnum(); return elseif (strcmpi(name,'CmMin')), enum=CmMinEnum(); return +elseif (strcmpi(name,'Gradient')), enum=GradientEnum(); return elseif (strcmpi(name,'Connectivity')), enum=ConnectivityEnum(); return elseif (strcmpi(name,'ControlParameter')), enum=ControlParameterEnum(); return Index: /issm/trunk/src/m/enum/TolXEnum.m =================================================================== --- /issm/trunk/src/m/enum/TolXEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/TolXEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=TolXEnum() -macro=244; +macro=245; Index: /issm/trunk/src/m/enum/VariableDescriptorsEnum.m =================================================================== --- /issm/trunk/src/m/enum/VariableDescriptorsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/VariableDescriptorsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=VariableDescriptorsEnum() -macro=245; +macro=246; Index: /issm/trunk/src/m/enum/VerboseEnum.m =================================================================== --- /issm/trunk/src/m/enum/VerboseEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/VerboseEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=VerboseEnum() -macro=246; +macro=247; Index: /issm/trunk/src/m/enum/WaitOnLockEnum.m =================================================================== --- /issm/trunk/src/m/enum/WaitOnLockEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/WaitOnLockEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=WaitOnLockEnum() -macro=247; +macro=248; Index: /issm/trunk/src/m/enum/YtsEnum.m =================================================================== --- /issm/trunk/src/m/enum/YtsEnum.m (revision 4046) +++ /issm/trunk/src/m/enum/YtsEnum.m (revision 4047) @@ -9,3 +9,3 @@ % macro=YtsEnum() -macro=248; +macro=249;