Changeset 27244

Timestamp:

08/26/22 21:25:22 (3 years ago)

Author:

Eric.Larour

Message:

CHG making sure we are not spreading a P0 thickness into P1 :)

File:

• issm/branches/trunk-larour-SLPS2022/src/c/analyses/MmemasstransportAnalysis.cpp (modified) (3 diffs)

issm/branches/trunk-larour-SLPS2022/src/c/analyses/MmemasstransportAnalysis.cpp

@@ -4 +4 @@
 #include "../classes/classes.h"
 #include "../classes/Inputs/TransientInput.h"
+#include "../classes/Inputs/TriaInput.h"
+#include "../classes/gauss/Gauss.h"
 #include "../shared/shared.h"
 #include "../modules/modules.h"
@@ -100 +102 @@
 void           MmemasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
-        /*retrieve thickness, ice and ocean levelsets from MmemasstransportAnalysis inputs in our element:*/
-
-        IssmDouble height,ice,ocean;
-        int       *doflist = NULL;
-
-        /*Fetch number of nodes and initialize values*/
-        int         numnodes = element->GetNumberOfNodes();
-        int         numdof   = numnodes*3;
-        IssmDouble* values   = xNew<IssmDouble>(numdof);
-
-        /*Get dof list and inputs */
-        element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-        Input* h_input=element->GetInput(MmemasstransportThicknessEnum); _assert_(h_input);
-        Input* i_input=element->GetInput(MmemasstransportMaskIceLevelsetEnum); _assert_(i_input);
-        Input* o_input=element->GetInput(MmemasstransportMaskOceanLevelsetEnum); _assert_(o_input);
-
-        /*Ok, we have the thickness in inputs, fill in solution */
-        Gauss* gauss=element->NewGauss();
-        for(int i=0;i<numnodes;i++){
-                gauss->GaussVertex(i);
-                h_input->GetInputValue(&height,gauss);
-                i_input->GetInputValue(&ice,gauss);
-                o_input->GetInputValue(&ocean,gauss);
-                values[3*i+0]=height;
-                values[3*i+1]=ice;
-                values[3*i+2]=ocean;
-        }
-
-        /*Add value to global vector*/
-        solution->SetValues(numdof,doflist,values,INS_VAL);
-
-        /*Free ressources:*/
-        delete gauss;
-        xDelete<int>(doflist);
-        xDelete<IssmDouble>(values);
+        /*do nothing:*/
+        return;
 }/*}}}*/
 void           MmemasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
@@ -141 +110 @@
 void           MmemasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
-        int         i,domaintype;
-        int*        doflist=NULL;
+        /*Fetch number of nodes and dof for this finite element*/
+        IssmDouble time;
+        IssmDouble ice[3];
+        IssmDouble ocean[3];
+        IssmDouble height;
+        int numnodes = element->GetNumberOfNodes();
+        
+        element->parameters->FindParam(&time,TimeEnum);
 
-        /*Fetch number of nodes and dof for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-        int numdof   = numnodes*3;
+        TriaInput* h_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportThicknessEnum,time)); _assert_(h_input);
+        TriaInput* i_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskIceLevelsetEnum,time)); _assert_(i_input);
+        TriaInput* o_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskOceanLevelsetEnum,time)); _assert_(o_input);
 
-        /*Fetch dof list and allocate solution vectors*/
-        element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
-        IssmDouble* values    = xNew<IssmDouble>(numdof);
-        IssmDouble* height        = xNew<IssmDouble>(numnodes);
-        IssmDouble* ice        = xNew<IssmDouble>(numnodes);
-        IssmDouble* ocean        = xNew<IssmDouble>(numnodes);
-
-        /*Use the dof list to index into the solution vector: */
-        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-
-        /*Retrieve h:*/
-        for(i=0;i<numnodes;i++){
-                height[i]=values[3*i+0];
-                ice[i]=values[3*i+1];
-                ocean[i]=values[3*i+2];
-
-                /*Check solution*/
-                if(xIsNan<IssmDouble>(height[i])) _error_("NaN found in ice thickness solution vector");
-                if(xIsInf<IssmDouble>(height[i])) _error_("Inf found in ice thickness solution vector");
-                if(xIsNan<IssmDouble>(ice[i])) _error_("NaN found in ice mask solution vector");
-                if(xIsInf<IssmDouble>(ice[i])) _error_("Inf found in ice mask solution vector");
-                if(xIsNan<IssmDouble>(ocean[i])) _error_("NaN found in ocean mask solution vector");
-                if(xIsInf<IssmDouble>(ocean[i])) _error_("Inf found in ocean mask solution vector");
-
+        Gauss* gauss=element->NewGauss();
+        for(int iv=0;iv<3;iv++){
+                gauss->GaussVertex(iv);
+                i_input->GetInputValue(&ice[iv],gauss);
+                o_input->GetInputValue(&ocean[iv],gauss);
         }
+        h_input->GetInputAverage(&height);
 
         /*Add thickness ice and ocean levelsets as inputs to the tria element: */
-        element->AddInput(ThicknessEnum,height,P0Enum); //very important, do not change the type of ThicknessEnum to P1 when it should be P0! 
+        element->AddInput(ThicknessEnum,&height,P0Enum); //very important, do not change the type of ThicknessEnum to P1 when it should be P0! 
         element->AddInput(MaskIceLevelsetEnum,ice,P1Enum);
         element->AddInput(MaskOceanLevelsetEnum,ocean,P1Enum);
-
-        /*Free ressources:*/
-        xDelete<IssmDouble>(height);
-        xDelete<IssmDouble>(ice);
-        xDelete<IssmDouble>(ocean);
-        xDelete<IssmDouble>(values);
-        xDelete<int>(doflist);
 
 }/*}}}*/