Changeset 16715

Timestamp:

11/12/13 12:09:40 (11 years ago)

Author:

seroussi

Message:

NEW: moved InputUpdateFromSolution in analyses for Stokes

Location:

issm/trunk-jpl/src/c

Files:

• analyses/L2ProjectionBaseAnalysis.cpp (modified) (1 diff)
• analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)
• analyses/StressbalanceAnalysis.h (modified) (1 diff)
• classes/Elements/Element.h (modified) (3 diffs)
• classes/Elements/Penta.cpp (modified) (2 diffs)
• classes/Elements/Penta.h (modified) (2 diffs)
• classes/Elements/Seg.h (modified) (3 diffs)
• classes/Elements/Tria.h (modified) (3 diffs)

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

@@ -59 +59 @@
 }/*}}}*/
 void L2ProjectionBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        int inputenum;
+
+        element->FindParam(&inputenum,InputToL2ProjectEnum);
+        element->InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
 }/*}}}*/

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

@@ -935 +935 @@
         switch(approximation){
                 case FSApproximationEnum: case NoneApproximationEnum:
-                        //InputUpdateFromSolutionFS(solution,element);
+                        InputUpdateFromSolutionFS(solution,element);
                         return;
                 case SSAApproximationEnum: case HOApproximationEnum: case SIAApproximationEnum:
@@ -1035 +1035 @@
 
 }/*}}}*/
+void StressbalanceAnalysis::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* values   = xNew<IssmDouble>(vnumdof+pnumdof);
+        IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
+        IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
+        IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
+        IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
+        IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
+
+        /*Prepare coordinate system list*/
+        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
+        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+        /*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++) values[i]        =solution[vdoflist[i]];
+        for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
+
+        /*Transform solution in Cartesian Space*/
+        element->TransformSolutionCoord(&values[0],&cs_list[0]);
+
+        /*Ok, we have vx and vy in values, fill in all arrays: */
+        for(i=0;i<vnumnodes;i++){
+                vx[i] = values[i*NDOF3+0];
+                vy[i] = values[i*NDOF3+1];
+                vz[i] = values[i*NDOF3+2];
+                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+                if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
+        }
+        for(i=0;i<pnumnodes;i++){
+                pressure[i] = values[vnumdof+i];
+                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
+        }
+
+        /*Recondition pressure and compute vel: */
+        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+        for(i = 0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
+        for(i = 0;i<vnumnodes;i++) vel[i]      = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+
+        /*Now, we have to move the previous inputs  to old 
+         * status, otherwise, we'll wipe them off: */
+        element->InputChangeName(VxEnum,VxPicardEnum);
+        element->InputChangeName(VyEnum,VyPicardEnum);
+        element->InputChangeName(VzEnum,VzPicardEnum);
+        element->InputChangeName(PressureEnum,PressurePicardEnum);
+
+        /*Add vx and vy as inputs to the tria element: */
+        element->AddInput(VxEnum,vx,P1Enum);
+        element->AddInput(VyEnum,vy,P1Enum);
+        element->AddInput(VzEnum,vz,P1Enum);
+        element->AddInput(VelEnum,vel,P1Enum);
+        element->AddInput(PressureEnum,pressure,P1Enum);
+
+        /*Free ressources:*/
+        xDelete<IssmDouble>(pressure);
+        xDelete<IssmDouble>(vel);
+        xDelete<IssmDouble>(vz);
+        xDelete<IssmDouble>(vy);
+        xDelete<IssmDouble>(vx);
+        xDelete<IssmDouble>(values);
+        xDelete<int>(vdoflist);
+        xDelete<int>(pdoflist);
+        xDelete<int>(cs_list);
+
+
+}/*}}}*/

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

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

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -50 +50 @@
                 virtual IssmDouble GetMaterialParameter(int enum_in)=0;
                 virtual void   TransformSolutionCoord(IssmDouble* values,int transformenum)=0;
+                virtual void   TransformSolutionCoord(IssmDouble* values,int* transformenum_list)=0;
                 virtual void    GetDofList(int** pdoflist,int approximation_enum,int setenum)=0;
                 virtual void    GetDofListVelocity(int** pdoflist,int setenum)=0;
@@ -56 +57 @@
                 virtual int    GetNodeIndex(Node* node)=0;
                 virtual int    GetNumberOfNodes(void)=0;
+                virtual int    GetNumberOfNodesVelocity(void)=0;
+                virtual int    GetNumberOfNodesPressure(void)=0;
                 virtual void   GetNodesSidList(int* sidlist)=0;
                 virtual void   GetNodesLidList(int* sidlist)=0;
@@ -85 +88 @@
                 virtual void   InputDuplicate(int original_enum,int new_enum)=0;
                 virtual void   InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code)=0;
+                virtual void   InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0;
 
                 virtual int    NodalValue(IssmDouble* pvalue, int index, int natureofdataenum)=0;

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -1359 +1359 @@
 }
 /*}}}*/
+/*FUNCTION Penta::GetNumberOfNodesPressure;{{{*/
+int Penta::GetNumberOfNodesPressure(void){
+        return this->NumberofNodesPressure();
+}
+/*}}}*/
+/*FUNCTION Penta::GetNumberOfNodesVelocity;{{{*/
+int Penta::GetNumberOfNodesVelocity(void){
+        return this->NumberofNodesVelocity();
+}
+/*}}}*/
 /*FUNCTION Penta::GetInput(int inputenum) {{{*/
 Input* Penta::GetInput(int inputenum){
@@ -3278 +3288 @@
 
         ::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum);
+
+}
+/*}}}*/
+/*FUNCTION Penta::TransformSolutionCoord{{{*/
+void Penta::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){
+
+        ::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum_list);
 
 }

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -91 +91 @@
                 void   GetNodesLidList(int* lidlist);
                 int    GetNumberOfNodes(void);
+                int    GetNumberOfNodesPressure(void);
+                int    GetNumberOfNodesVelocity(void);
                 IssmDouble GetMaterialParameter(int enum_in);
                 void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
@@ -242 +244 @@
                 void             SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
                 void           TransformSolutionCoord(IssmDouble* values,int transformenum);
+                void           TransformSolutionCoord(IssmDouble* values,int* transformenum_list);
 
                 #ifdef _HAVE_STRESSBALANCE_

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -54 +54 @@
                 /*}}}*/
                 /*Update virtual functions resolution: {{{*/
-                void  InputUpdateFromSolution(IssmDouble* solutiong){_error_("not implemented yet");};
+                void  InputUpdateFromSolution(IssmDouble* solution){_error_("not implemented yet");};
+                void  InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
                 void  InputUpdateFromVector(IssmDouble* vector, int name, int type){_error_("not implemented yet");};
 #ifdef _HAVE_DAKOTA_
@@ -93 +94 @@
                 void        GetNodesLidList(int* lidlist){_error_("not implemented yet");};
                 int         GetNumberOfNodes(void){_error_("not implemented yet");};
+                int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
+                int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
                 void        GetVerticesCoordinates(IssmDouble** pxyz_list){_error_("not implemented yet");};
                 int         Sid(){_error_("not implemented yet");};
@@ -143 +146 @@
                 IssmDouble  TimeAdapt(){_error_("not implemented yet");};
                 void        TransformSolutionCoord(IssmDouble* values,int transformenum){_error_("not implemented yet");};
+                void        TransformSolutionCoord(IssmDouble* values,int* transformenum_list){_error_("not implemented yet");};
                 void UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};
                 void UpdateConstraintsExtrudeFromTop(){_error_("not implemented");};

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -89 +89 @@
                 void        GetNodesLidList(int* lidlist);
                 int         GetNumberOfNodes(void);
+                int         GetNumberOfNodesPressure(void){_error_("not implemented yet");};
+                int         GetNumberOfNodesVelocity(void){_error_("not implemented yet");};
                 int         Sid();
                 bool        IsOnBed();
@@ -257 +259 @@
                 void             InputChangeName(int enum_type,int enum_type_old);
                 void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
+                void             InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
                 void             InputUpdateFromSolutionMasstransport(IssmDouble* solution);
                 bool             IsInput(int name);
@@ -264 +267 @@
                 void             SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
                 void           TransformSolutionCoord(IssmDouble* values,int transformenum);
+                void           TransformSolutionCoord(IssmDouble* values,int* transformenum_list){_error_("not implemented yet");};
                 #ifdef _HAVE_STRESSBALANCE_
                 ElementMatrix* CreateKMatrixStressbalanceSSA(void);