Changeset 16752

Timestamp:

11/13/13 21:16:11 (11 years ago)

Author:

seroussi

Message:

CHG: added InputUpdateFromSolution for Adjoint Stokes

Location:

issm/trunk-jpl/src/c/analyses

Files:

• AdjointHorizAnalysis.cpp (modified) (2 diffs)
• AdjointHorizAnalysis.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

@@ -30 +30 @@
 }/*}}}*/
 void AdjointHorizAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+        int approximation;
+        element->GetInputValue(&approximation,ApproximationEnum);
+        if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
+                InputUpdateFromSolutionFS(solution,element);
+        }
+        else{
+                InputUpdateFromSolutionHoriz(solution,element);
+        }
+}/*}}}*/
+void AdjointHorizAnalysis::InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element){/*{{{*/
         int  i;
         int* doflist=NULL;
@@ -69 +79 @@
         xDelete<int>(doflist);
 }/*}}}*/
+void AdjointHorizAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
+        int          i;
+        int*         vdoflist=NULL;
+        int*         pdoflist=NULL;
+        IssmDouble   FSreconditioning;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int vnumnodes = element->GetNumberOfNodesVelocity();
+        int pnumnodes = element->GetNumberOfNodesPressure();
+        int vnumdof   = vnumnodes*3;
+        int pnumdof   = pnumnodes*1;
+
+        /*Initialize values*/
+        IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
+        IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
+        IssmDouble* lambdax = xNew<IssmDouble>(vnumnodes);
+        IssmDouble* lambday = xNew<IssmDouble>(vnumnodes);
+        IssmDouble* lambdaz = xNew<IssmDouble>(vnumnodes);
+        IssmDouble* lambdap = xNew<IssmDouble>(pnumnodes);
+
+        /*Get dof list: */
+        element->GetDofListVelocity(&vdoflist,GsetEnum);
+        element->GetDofListPressure(&pdoflist,GsetEnum);
+
+        /*Use the dof list to index into the solution vector: */
+        for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
+        for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
+
+        /*Transform solution in Cartesian Space*/
+        element->TransformSolutionCoord(&vvalues[0],XYZEnum);
+
+        /*fill in all arrays: */
+        for(i=0;i<vnumnodes;i++){
+                lambdax[i] = vvalues[i*NDOF3+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
+                lambday[i] = vvalues[i*NDOF3+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
+                lambdaz[i] = vvalues[i*NDOF3+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
+        }
+        for(i=0;i<pnumnodes;i++){
+                lambdap[i] = pvalues[i]; if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
+        }
+
+        /*Recondition pressure and compute vel: */
+        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+        for(i=0;i<pnumnodes;i++) lambdap[i]=lambdap[i]*FSreconditioning;
+
+        /*Add vx and vy as inputs to the tria element: */
+        element->AddInput(AdjointxEnum,lambdax,P1Enum);
+        element->AddInput(AdjointyEnum,lambday,P1Enum);
+        element->AddInput(AdjointzEnum,lambdaz,P1Enum);
+        element->AddInput(AdjointpEnum,lambdap,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<int>(vdoflist);
+        xDelete<int>(pdoflist);
+        xDelete<IssmDouble>(lambdap);
+        xDelete<IssmDouble>(lambdaz);
+        xDelete<IssmDouble>(lambday);
+        xDelete<IssmDouble>(lambdax);
+        xDelete<IssmDouble>(vvalues);
+        xDelete<IssmDouble>(pvalues);
+}/*}}}*/

issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h

@@ -20 +20 @@
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
+                void InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element);
+                void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
 };
 #endif