Changeset 4651

Timestamp:

07/19/10 11:24:03 (15 years ago)

Author:

Mathieu Morlighem

Message:

Fixed MacAyeal-Hutter and some debugging in CM Stokes. Also removed a couple of useless files

Location:

issm/trunk/src

Files:

• c/objects/Elements/Tria.cpp (modified) (2 diffs)
• c/objects/Elements/Tria.h (modified) (1 diff)
• c/solutions/diagnostic_core.cpp (modified) (1 diff)
• c/solutions/stokescontrolinit.cpp (modified) (1 diff)
• m/solutions/ControlInitialization.m (deleted)
• m/solutions/diagnostic_core.m (modified) (1 diff)
• m/solutions/gradjcompute_core.m (deleted)
• m/solutions/stokescontrolinit.m (modified) (1 diff)

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

@@ -1034 +1034 @@
         /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
         if (analysis_type==DiagnosticHorizAnalysisEnum)
-                GetSolutionFromInputsDiagnosticHoriz(solution);
+         GetSolutionFromInputsDiagnosticHoriz(solution);
         else if (analysis_type==AdjointHorizAnalysisEnum)
          GetSolutionFromInputsAdjointHoriz(solution);
+        else if (analysis_type==DiagnosticHutterAnalysisEnum)
+         GetSolutionFromInputsDiagnosticHutter(solution);
         else
          ISSMERROR("analysis: %s not supported yet",EnumAsString(analysis_type));
@@ -5352 +5354 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
+void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
+
+        int i;
+
+        const int    numvertices=3;
+        const int    numdofpervertex=2;
+        const int    numdof=numdofpervertex*numvertices;
+        double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
+
+        int          doflist[numdof];
+        double       values[numdof];
+        double       vx;
+        double       vy;
+
+        int          dummy;
+
+        /*Get dof list: */
+        GetDofList(&doflist[0],&dummy);
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        /*P1 element only for now*/
+        for(i=0;i<numvertices;i++){
+
+                /*Recover vx and vy*/
+                inputs->GetParameterValue(&vx,&gauss[i][0],VxEnum);
+                inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum);
+                values[i*numdofpervertex+0]=vx;
+                values[i*numdofpervertex+1]=vy;
+        }
+
+        /*Add value to global vector*/
+        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
+
+}
+/*}}}*/
 /*FUNCTION Tria::GradjDragStokes {{{1*/
 void  Tria::GradjDragStokes(Vec gradient){

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

@@ -154 +154 @@
                 void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
                 void      GetSolutionFromInputsAdjointHoriz(Vec solution);
+                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);
                 void      GradjDragStokes(Vec gradient);
                 void      InputUpdateFromSolutionAdjointHoriz( double* solution);

issm/trunk/src/c/solutions/diagnostic_core.cpp

@@ -51 +51 @@
                         
                 if(verbose)_printf_("%s\n"," computing hutter velocities...");
+
+                //Take the last velocity into account so that the velocity on the MacAyeal domain is not zero
+                if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHutterAnalysisEnum);
+
                 femmodel->SetCurrentConfiguration(DiagnosticHutterAnalysisEnum);
                 solver_linear(NULL,femmodel);

issm/trunk/src/c/solutions/stokescontrolinit.cpp

@@ -42 +42 @@
         femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
         solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
-
-        //vertical velocity
+        femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
         solver_linear(NULL,femmodel);
 

issm/trunk/src/m/solutions/diagnostic_core.m

@@ -35 +35 @@
 
                 displaystring(verbose,'\n%s',['computing hutter velocities...']);
+
+                %Take the last velocity into account so that the velocity on the MacAyeal domain is not zero
+                if(ismacayealpattyn)femmodel=ResetBoundaryConditions(femmodel,DiagnosticHutterAnalysisEnum); end 
+
                 femmodel=SetCurrentConfiguration(femmodel,DiagnosticHutterAnalysisEnum);
                 femmodel=solver_linear(femmodel);

issm/trunk/src/m/solutions/stokescontrolinit.m

@@ -21 +21 @@
         femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
         femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads);
-
-        %vertical velocity
+        femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
         femmodel=solver_linear(femmodel);