Changeset 16908

Timestamp:

11/23/13 20:22:24 (11 years ago)

Author:

seroussi

Message:

BUG: fixed SSAHO PVector Coupling Viscous

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (5 diffs)
• classes/Elements/Penta.cpp (modified) (1 diff)

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

@@ -3096 +3096 @@
 
         /*Constants*/
-        int numnodes = 2*element->GetNumberOfNodes();
+        int numnodes    = element->GetNumberOfNodes();
         int numdofm     = 1 *numnodes; //*2/2
         int numdofp     = 2 *numnodes;
@@ -3111 +3111 @@
         IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
         IssmDouble* Ke_gg_gaussian = xNew<IssmDouble>(numdofp*numdofm);
-        Node        *node_list[numnodes];
-        int         cs_list[numnodes];
+        Node        *node_list[2*numnodes];
+        int*         cs_list= xNew<int>(2*numnodes);
 
         /*Find penta on bed as HO must be coupled to the dofs on the bed: */
@@ -3150 +3150 @@
                 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
                 this->GetBSSAHO(B, element,xyz_list, gauss);
-                //basaltria->GetBprimeSSA(Bprime, xyz_list, gauss_tria); /FIXME
                 this->GetBSSAprime(Bprime, basaltria,xyz_list, gauss_tria); 
                 element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
@@ -3164 +3163 @@
                                         Ke_gg_gaussian,0);
 
-                for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i*numdofp+j]+=Ke_gg_gaussian[i*numdofp+j];
+                for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i*numdofm+j]+=Ke_gg_gaussian[i*numdofm+j];
         } 
-        for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofp+j];
-        for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofp+i];
+        for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
+        for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
 
         /*Transform Coordinate System*/
-        element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+        element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
 
         /*Clean-up and return*/
         basaltria->DeleteMaterials(); delete basaltria;
+        
         delete gauss;
         delete gauss_tria;
+        xDelete<IssmDouble>(Ke_gg);
+        xDelete<IssmDouble>(Ke_gg_gaussian);
+        xDelete<int>(cs_list);
         return Ke;
 
@@ -3211 +3214 @@
         /*Build B: */
         for(int i=0;i<numnodes;i++){
-                B[2*numnodes*0+2*i+0] = dbasis[0*3+i];
+                B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
                 B[2*numnodes*0+2*i+1] = 0.;
                 B[2*numnodes*1+2*i+0] = 0.;
-                B[2*numnodes*1+2*i+1] = dbasis[1*3+i];
-                B[2*numnodes*2+2*i+0] = .5*dbasis[1*3+i];
-                B[2*numnodes*2+2*i+1] = .5*dbasis[0*3+i];
+                B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
+                B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
+                B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
         }
 

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

@@ -1328 +1328 @@
 /*FUNCTION Penta::GetNode(int node_number) {{{*/
 Node* Penta::GetNode(int node_number){
+        _assert_(node_number>=0); 
+        _assert_(node_number<this->NumberofNodes()); 
         return this->nodes[node_number];
 }