Changeset 4353

Timestamp:

06/30/10 14:17:09 (15 years ago)

Author:

Mathieu Morlighem

Message:

fixing control methods

Location:

issm/trunk/src/c

Files:

• objects/Elements/Penta.cpp (modified) (3 diffs)
• objects/Elements/Tria.cpp (modified) (7 diffs)
• objects/Elements/Tria.h (modified) (1 diff)
• objects/Inputs/IntInput.cpp (modified) (1 diff)
• shared/Exceptions/Exceptions.cpp (modified) (1 diff)
• solutions/adjoint_core.cpp (modified) (1 diff)
• solutions/control_core.cpp (modified) (4 diffs)
• solutions/gradient_core.cpp (modified) (1 diff)
• solutions/objectivefunctionC.cpp (modified) (2 diffs)

issm/trunk/src/c/objects/Elements/Penta.cpp ¶

@@ -290 +290 @@
 /*FUNCTION Penta::InputUpdateFromConstant(bool value, int name);{{{1*/
 void  Penta::InputUpdateFromConstant(bool constant, int name){
-        /*Nothing updated for now*/
+
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new BoolInput(name,constant));
 }
 /*}}}*/
 /*FUNCTION Penta::InputUpdateFromConstant(double value, int name);{{{1*/
 void  Penta::InputUpdateFromConstant(double constant, int name){
-        /*Nothing updated for now*/
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new DoubleInput(name,constant));
 }
 /*}}}*/
 /*FUNCTION Penta::InputUpdateFromConstant(int value, int name);{{{1*/
 void  Penta::InputUpdateFromConstant(int constant, int name){
-        /*Nothing updated for now*/
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new IntInput(name,constant));
 }
 /*}}}*/
@@ -2247 +2260 @@
 /*FUNCTION Penta::IsInput{{{1*/
 bool Penta::IsInput(int name){
-        if (name==ThicknessEnum ||
+        if (
+                                name==ThicknessEnum ||
                                 name==SurfaceEnum ||
                                 name==BedEnum ||
@@ -2262 +2276 @@
                                 name==TemperatureAverageEnum ||
                                 name==RheologyBEnum ||
-                                name==RheologyNEnum) {
+                                name==RheologyNEnum ||
+                                name==FitEnum
+                                ) {
                 return true;
         }

issm/trunk/src/c/objects/Elements/Tria.cpp ¶

@@ -302 +302 @@
 /*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{1*/
 void  Tria::InputUpdateFromConstant(int constant, int name){
-        /*Nothing updated for now*/
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new IntInput(name,constant));
 }
 /*}}}*/
 /*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{1*/
 void  Tria::InputUpdateFromConstant(double constant, int name){
-        /*Nothing updated for now*/
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new DoubleInput(name,constant));
 }
 /*}}}*/
 /*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{1*/
 void  Tria::InputUpdateFromConstant(bool constant, int name){
-        /*Nothing updated for now*/
+        /*Check that name is an element input*/
+        if (!IsInput(name)) return;
+
+        /*update input*/
+        this->inputs->AddInput(new BoolInput(name,constant));
 }
 /*}}}*/
@@ -329 +341 @@
         else if (analysis_type==DiagnosticHorizAnalysisEnum){
                 InputUpdateFromSolutionDiagnosticHoriz( solution);
+        }
+        else if (analysis_type==AdjointAnalysisEnum){
+                InputUpdateFromSolutionAdjoint( solution);
         }
         else if (analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum){
@@ -721 +736 @@
         double scale=0;
         double S=0;
-        double fit=-1;
+        int    fit=-1;
 
         /*retrieve some parameters: */
@@ -1887 +1902 @@
         double scaley=1;
         double S=0;
-        double fit=-1;
+        int    fit=-1;
 
         /*inputs: */
@@ -2120 +2135 @@
         /*inputs: */
         bool onwater;
-        int  fit;
+        int fit;
 
         /*retrieve inputs :*/
@@ -5319 +5334 @@
 }
 /*}}}*/
+/*FUNCTION Tria::InputUpdateFromSolutionAdjoint {{{1*/
+void  Tria::InputUpdateFromSolutionAdjoint(double* solution){
+
+        int i;
+
+        const int    numvertices=3;
+        const int    numdofpervertex=2;
+        const int    numdof=numdofpervertex*numvertices;
+
+        int          doflist[numdof];
+        double       values[numdof];
+        double       lambdax[numvertices];
+        double       lambday[numvertices];
+
+        int          dummy;
+
+        /*Get dof list: */
+        GetDofList(&doflist[0],&dummy);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<numdof;i++){
+                values[i]=solution[doflist[i]];
+        }
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        for(i=0;i<numvertices;i++){
+                lambdax[i]=values[i*numdofpervertex+0];
+                lambday[i]=values[i*numdofpervertex+1];
+        }
+
+        /*Add vx and vy as inputs to the tria element: */
+        this->inputs->AddInput(new TriaVertexInput(AdjointxEnum,lambdax));
+        this->inputs->AddInput(new TriaVertexInput(AdjointyEnum,lambday));
+
+}
+/*}}}*/
 /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
 void  Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
@@ -5473 +5524 @@
 /*FUNCTION Tria::IsInput{{{1*/
 bool Tria::IsInput(int name){
-        if (name==SurfaceSlopeXEnum ||
-                                name==SurfaceSlopeYEnum){
+        if (
+                                name==SurfaceSlopeXEnum ||
+                                name==SurfaceSlopeYEnum ||
+                                name==FitEnum
+                ){
                 return true;
         }

issm/trunk/src/c/objects/Elements/Tria.h ¶

@@ -152 +152 @@
                 void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
                 void      GradjDragStokes(Vec gradient);
+                void      InputUpdateFromSolutionAdjoint( double* solution);
                 void      InputUpdateFromSolutionDiagnosticHoriz( double* solution);
                 void      InputUpdateFromSolutionSlopeCompute( double* solution);

issm/trunk/src/c/objects/Inputs/IntInput.cpp ¶

@@ -191 +191 @@
 /*}}}*/
 /*FUNCTION IntInput::GetParameterValue(double* pvalue){{{1*/
-void IntInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");}
+void IntInput::GetParameterValue(double* pvalue){
+        ISSMERROR("IntInput cannot return a double in parallel");
+}
 /*}}}*/
 /*FUNCTION IntInput::GetParameterValue(double* pvalue,Node* node){{{1*/

issm/trunk/src/c/shared/Exceptions/Exceptions.cpp ¶

@@ -42 +42 @@
         else{
                 #ifdef _PARALLEL_
-                        printf("\n[%i] %s%s%s%i\n",my_rank,"??? Error using ==> ",file_name.c_str()," at ",file_line);
-                        printf("[%i] %s%s%s\n\n",my_rank,function_name.c_str()," error message: ",what());
+                        printf("\n[%i] ??? Error using ==> %s:%i\n",my_rank,file_name.c_str(),file_line);
+                        printf("[%i] %s error message: %s\n\n",my_rank,function_name.c_str(),what());
                 #else
-                        printf("\n%s%s%s%i\n","??? Error using ==> ",file_name.c_str()," at ",file_line);
-                        printf("%s%s%s\n\n",function_name.c_str()," error message: ",what());
+                        printf("\n??? Error using ==> %s at %i\n",file_name.c_str(),file_line);
+                        printf("%s error message: %s\n\n",function_name.c_str(),what());
                 #endif
         }

issm/trunk/src/c/solutions/adjoint_core.cpp ¶

@@ -11 +11 @@
 #include "../include/include.h"
 #include "../solvers/solvers.h"
-
 
 void adjoint_core(FemModel* femmodel){

issm/trunk/src/c/solutions/control_core.cpp ¶

@@ -19 +19 @@
         int     verbose=0;
         int     control_type;
-        int     control_steady;
+        bool    control_steady;
         int     nsteps;
         double  eps_cm;
@@ -25 +25 @@
         double  cm_min;
         double  cm_max;
-        int     cm_gradient;
+        bool    cm_gradient;
         int     dim;
 
@@ -48 +48 @@
         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(&fit,NULL,FitEnum);
+        femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
+        femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum);
+        femmodel->parameters->FindParam(&cm_jump,NULL,CmJumpEnum);
         femmodel->parameters->FindParam(&eps_cm,EpsCmEnum);
         femmodel->parameters->FindParam(&tolx,TolXEnum);
@@ -72 +72 @@
 
                 _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
-                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,fit[n],FitEnum);
+                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)fit[n],FitEnum);
                 
                 /*In case we are running a steady state control method, compute new temperature field using new parameter * distribution: */

issm/trunk/src/c/solutions/gradient_core.cpp ¶

@@ -16 +16 @@
         
         /*parameters: */
-        int control_steady;
+        bool control_steady;
         int control_type;
         int verbose;

issm/trunk/src/c/solutions/objectivefunctionC.cpp ¶

@@ -36 +36 @@
         int     control_type;
         int     analysis_type;
-        int     control_steady;
+        bool    control_steady;
         bool    conserve_loads=true;
         
@@ -45 +45 @@
         n=optargs->n;
 
-        femmodel->parameters->FindParam(&optscal,NULL,NULL,OptScalEnum);
-        femmodel->parameters->FindParam(&fit,NULL,NULL,FitEnum);
+        femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum);
+        femmodel->parameters->FindParam(&fit,NULL,FitEnum);
         femmodel->parameters->FindParam(&cm_min,CmMinEnum);
         femmodel->parameters->FindParam(&cm_max,CmMaxEnum);