Changeset 15872

Timestamp:

08/22/13 08:31:52 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing some stuff for Johannes

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (5 diffs)
• PentaRef.cpp (modified) (1 diff)
• PentaRef.h (modified) (1 diff)

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

@@ -5026 +5026 @@
 
         /*Fetch indices of basal nodes for this finite element*/
-        BasalNodeIndices(&numindices,&indices,this->VelocityInterpolation());
+        BasalNodeIndices(&numindices,&indices,this->element_type);
 
         /*Get parameters and inputs: */
@@ -5062 +5062 @@
         /* Intermediaries */
         bool        isenthalpy;
-        int         i,j,analysis_type;
+        int         i,j,analysis_type,numindices,numindicesup;
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -5077 +5077 @@
         IssmDouble  vx,vy,vz;
         IssmDouble  dt;
-        Friction   *friction                 = NULL;
-        GaussPenta *gauss                    = NULL;
-        GaussPenta *gaussup                  = NULL;
+        int        *indices   = NULL;
+        int        *indicesup = NULL;
+        Friction   *friction  = NULL;
 
         /* Only compute melt rates at the base of grounded ice*/
@@ -5104 +5104 @@
         friction=new Friction("3d",inputs,matpar,analysis_type);
 
-        /* Start  looping on the number of gaussian points: */
-        gauss=new GaussPenta(0,1,2,2);
-        gaussup=new GaussPenta(3,4,5,2);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
-                gauss->GaussPoint(ig);
-                gaussup->GaussPoint(ig);
+        /*Fetch indices of basal and surface nodes for this finite element*/
+        BasalNodeIndices(&numindices,&indices,this->element_type);
+        SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
+        _assert_(numindices==numindicesup); 
+
+        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+        GaussPenta* gauss=new GaussPenta();
+        GaussPenta* gaussup=new GaussPenta();
+        for(i=0;i<numindices;i++){
+                gauss->GaussNode(this->element_type,indices[i]);
+                gaussup->GaussNode(this->element_type,indicesup[i]);
 
                 watercolumn_input->GetInputValue(&watercolumn, gauss);
@@ -5161 +5166 @@
                 // TODO: feed meltrate & watercolumn back to model
         }  
+
+        /*Clean up and return*/
+        xDelete<int>(indices);
+        xDelete<int>(indicesup);
+        delete gauss;
+        delete gaussup;
 }
 /*}}}*/

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

@@ -2118 +2118 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::SurfaceNodeIndices{{{*/
+void PentaRef::SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement){
+
+        /*Output*/
+        int  numindices;
+        int* indices = NULL;
+
+        switch(finiteelement){
+                case P1Enum: case P1DGEnum:
+                        numindices = 3;
+                        indices    = xNew<int>(numindices);
+                        indices[0] = 3;
+                        indices[1] = 4;
+                        indices[2] = 5;
+                        break;
+                case P1bubbleEnum: case P1bubblecondensedEnum:
+                        numindices = 3;
+                        indices    = xNew<int>(numindices);
+                        indices[0] = 3;
+                        indices[1] = 4;
+                        indices[2] = 5;
+                        break;
+                case P2xP1Enum:
+                        numindices = 6;
+                        indices    = xNew<int>(numindices);
+                        indices[0] = 3;
+                        indices[1] = 4;
+                        indices[2] = 5;
+                        indices[3] = 9;
+                        indices[4] = 10;
+                        indices[5] = 11;
+                        break;
+                case P1xP2Enum:
+                        numindices = 3;
+                        indices    = xNew<int>(numindices);
+                        indices[0] = 3;
+                        indices[1] = 4;
+                        indices[2] = 5;
+                        return;
+                case P2Enum:
+                        numindices = 6;
+                        indices    = xNew<int>(numindices);
+                        indices[0] = 3;
+                        indices[1] = 4;
+                        indices[2] = 5;
+                        indices[3] = 12;
+                        indices[4] = 13;
+                        indices[5] = 14;
+                        break;
+                default:
+                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+        }
+
+        /*Assign output pointer*/
+        *pnumindices = numindices;
+        *pindices    = indices;
+}
+/*}}}*/

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

@@ -70 +70 @@
 
                 void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
+                void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
                 int  NumberofNodes(void);
                 int  NumberofNodesVelocity(void);