Changeset 15770

Timestamp:

08/09/13 10:44:43 (12 years ago)

Author:

seroussi

Message:

BUG: starting to debug tiling HO/FS for ice sheets

File:

• issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (modified) (5 diffs)

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

@@ -6990 +6990 @@
 
         /*Constants*/
-        const int numnodes  = 2 *NUMVERTICES;
         const int numdof    = NUMVERTICES *NDOF4;
         const int numdofm   = NUMVERTICES *NDOF2;
@@ -7016 +7015 @@
         Friction*  friction=NULL;
         GaussPenta *gauss=NULL;
-        Node       *node_list[numnodes];
-        int         cs_list[numnodes];
 
         /*If on water or not FS, skip stiffness: */
         inputs->GetInputValue(&approximation,ApproximationEnum);
         if(IsFloating() || !IsOnBed()) return NULL;
+
+        int vnumnodes = this->NumberofNodesVelocity();
+        int pnumnodes = this->NumberofNodesPressure();
+        int numnodes  = 2*vnumnodes-1+pnumnodes;
+
+        Node       *node_list[numnodes];
+        int         cs_list[numnodes];
+
         ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum);
-        ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
+        ElementMatrix* Ke2=new ElementMatrix(this->nodes,numnodes,this->parameters,FSvelocityEnum);
         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
         delete Ke1; delete Ke2;
 
         /*Prepare node list*/
-        for(i=0;i<NUMVERTICES;i++){
+        for(i=0;i<numnodes+NUMVERTICES;i++){
                 node_list[i+0*NUMVERTICES] = this->nodes[i];
                 node_list[i+1*NUMVERTICES] = this->nodes[i];
@@ -7097 +7102 @@
 
         /*Clean up and return*/
+        xDelete<int>(cs_list);
         delete gauss;
         delete friction;
@@ -8352 +8358 @@
         inputs->GetInputValue(&approximation,ApproximationEnum);
         if(approximation!=HOFSApproximationEnum) return NULL;
-        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
+
+        int vnumnodes = this->NumberofNodesVelocity();
+        int pnumnodes = this->NumberofNodesPressure();
+        int numnodes  = vnumnodes+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;
+
+        ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,FSvelocityEnum);
 
         /*Retrieve all inputs and parameters*/
@@ -8386 +8402 @@
 
                 for(i=0;i<NUMVERTICES2D;i++){
-                        pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
-                        pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
-                        pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
+                        pe->values[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
+                        pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
+                        pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
                 }
         }
 
         /*Transform coordinate system*/
-        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
+        TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
         /*Clean up and return*/
+        xDelete<int>(cs_list);
         delete gauss;
         delete friction;