Changeset 16756

Timestamp:

11/14/13 09:43:52 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: extended InputUpdateFromSolution for FS 2d

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (6 diffs)
• classes/Elements/Tria.cpp (modified) (1 diff)
• classes/Elements/Tria.h (modified) (1 diff)
• modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp (modified) (1 diff)

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

@@ -1130 +1130 @@
 void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
 
-        int          i;
+        int          i,dim,meshtype;
         int*         vdoflist=NULL;
         int*         pdoflist=NULL;
         IssmDouble   FSreconditioning;
 
+        element->FindParam(&meshtype,MeshTypeEnum);
+        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+        switch(meshtype){
+                case Mesh2DverticalEnum: dim = 2; break;
+                case Mesh3DEnum:         dim = 3; break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
         /*Fetch number of nodes and dof for this finite element*/
         int vnumnodes = element->GetNumberOfNodesVelocity();
         int pnumnodes = element->GetNumberOfNodesPressure();
-        int vnumdof   = vnumnodes*3;
+        int vnumdof   = vnumnodes*dim;
         int pnumdof   = pnumnodes*1;
 
@@ -1151 +1159 @@
         /*Prepare coordinate system list*/
         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
-        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
+        if(dim==2){
+                for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+        }
+        else{
+                for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
+        }
         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
 
@@ -1167 +1180 @@
         /*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];
+                vx[i] = values[i*dim+0];
+                vy[i] = values[i*dim+1];
                 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");
+
+                if(dim==3){
+                        vz[i] = values[i*dim+2];
+                        if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
+                }
         }
         for(i=0;i<pnumnodes;i++){
@@ -1180 +1196 @@
 
         /*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]);
+        for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
+        if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+        else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
 
         /*Now, we have to move the previous inputs  to old 
@@ -1188 +1204 @@
         element->InputChangeName(VxEnum,VxPicardEnum);
         element->InputChangeName(VyEnum,VyPicardEnum);
-        element->InputChangeName(VzEnum,VzPicardEnum);
         element->InputChangeName(PressureEnum,PressurePicardEnum);
+        if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
 
         /*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);
+        if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
 
         /*Free ressources:*/
@@ -1208 +1224 @@
         xDelete<int>(pdoflist);
         xDelete<int>(cs_list);
-
-
 }/*}}}*/

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

@@ -2847 +2847 @@
 
         ::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum);
+
+}
+/*}}}*/
+/*FUNCTION Tria::TransformSolutionCoord{{{*/
+void Tria::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){
+
+        ::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum_list);
 
 }

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

@@ -271 +271 @@
                 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");};
+                void           TransformSolutionCoord(IssmDouble* values,int* transformenum_list);
                 #ifdef _HAVE_STRESSBALANCE_
                 ElementMatrix* CreateKMatrixStressbalanceSSA(void);

issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp

@@ -28 +28 @@
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                //analysis->InputUpdateFromSolution(solution,element);
-                element->InputUpdateFromSolution(solution);
+                analysis->InputUpdateFromSolution(solution,element);
+                //element->InputUpdateFromSolution(solution);
         }
         delete analysis;