Changeset 17315

Timestamp:

02/20/14 08:55:05 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on SSA 1d

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (5 diffs)
• classes/Elements/Seg.cpp (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/gauss/GaussSeg.cpp (modified) (1 diff)

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

@@ -16 +16 @@
 
         switch(approximation){
-                case SSAApproximationEnum:  numdofs =2; break;
+                case SSAApproximationEnum:
+                         switch(meshtype){
+                                 case Mesh3DEnum:         numdofs=2; break;
+                                 case Mesh2DverticalEnum: numdofs=1; break;
+                                 default: _error_("mesh type not supported yet");
+                         }
+                         break;
                 case L1L2ApproximationEnum: numdofs =2; break;
                 case HOApproximationEnum:   
@@ -1422 +1428 @@
 
         /*Transform coordinate system*/
-        element->TransformLoadVectorCoord(pe,XYEnum);
+        if(dim==2) element->TransformLoadVectorCoord(pe,XYEnum);
 
         /*Clean up and return*/
@@ -1496 +1502 @@
 
         /*Transform coordinate system*/
-        element->TransformLoadVectorCoord(pe,XYEnum);
+        if(dim==2) element->TransformLoadVectorCoord(pe,XYEnum);
 
         /*Clean up and return*/
@@ -1696 +1702 @@
         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
         for(i=0;i<numnodes;i++){
-                vx[i]=values[i*2+0];
+                vx[i]=values[i*dim+0];
                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
 
                 if(dim==2){
-                        vy[i]=values[i*2+1];
+                        vy[i]=values[i*dim+1];
                         if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
                 }
@@ -1711 +1717 @@
         }
         else{
-                element->GetInputListOnNodes(&vy[0],VyEnum,0.);
+                basalelement->GetInputListOnNodes(&vy[0],VyEnum,0.);
                 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
         }

issm/trunk-jpl/src/c/classes/Elements/Seg.cpp ¶

@@ -227 +227 @@
 }
 /*}}}*/
+/*FUNCTION Seg::JacobianDeterminantSurface{{{*/
+void Seg::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
+
+        *pJdet = 1.;
+
+}
+/*}}}*/
 /*FUNCTION Seg::NewGauss(){{{*/
 Gauss* Seg::NewGauss(void){

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

@@ -96 +96 @@
                 void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
                 void        JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
-                void        JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
+                void        JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
                 void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
                 void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};

issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp ¶

@@ -133 +133 @@
 void GaussSeg::GaussNode(int finiteelement,int iv){
 
-        _error_("not supported");
+        /*in debugging mode: check that the default constructor has been called*/
+        _assert_(numgauss==-1);
+
+        /*update static arrays*/
+        switch(finiteelement){
+                case P1Enum: case P1DGEnum:
+                        switch(iv){
+                                case 0: coord1=-1.; break;
+                                case 1: coord1=+1.; break;
+                                default: _error_("node index should be in [0 1]");
+                        }
+                        break;
+                default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
+        }
 
 }