source: issm/oecreview/Archive/16554-17801/ISSM-16975-16976.diff

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

@@ -3214 +3214 @@
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/
 
         /*compute all stiffness matrices for this element*/
+        printf("-------------- file: StressbalanceAnalysis.cpp line: %i\n",__LINE__); 
         ElementMatrix* Ke1=CreateKMatrixFS(element);
+        printf("-------------- file: StressbalanceAnalysis.cpp line: %i\n",__LINE__); 
         int indices[3]={18,19,20};
         Ke1->StaticCondensation(3,&indices[0]);
+        printf("-------------- file: StressbalanceAnalysis.cpp line: %i\n",__LINE__); 
         int init = element->FiniteElement();
         element->SetTemporaryElementType(P1Enum); 
+        printf("-------------- file: StressbalanceAnalysis.cpp line: %i\n",__LINE__); 
         ElementMatrix* Ke2=CreateKMatrixSSA3d(element);
+        printf("-------------- file: StressbalanceAnalysis.cpp line: %i\n",__LINE__); 
         element->SetTemporaryElementType(init); 
         ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element);
+        printf("-------------- file: StressbalanceAnalysis.cpp line: %i\n",__LINE__); 
         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
         /*clean-up and return*/
@@ -3535 +3541 @@
                 cs_list[2*vnumnodes-1+i]   = PressureEnum;
         }
 
-        ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
+        ElementMatrix* Ke1=element->NewElementMatrixCoupling(6,SSAApproximationEnum);
         ElementMatrix* Ke2=element->NewElementMatrix(FSvelocityEnum);
         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
         delete Ke1; delete Ke2;
@@ -3641 +3647 @@
         int* cs_list     = xNew<int>(2*vnumnodes-1+pnumnodes);
         Node **node_list = xNew<Node*>(2*vnumnodes-1+pnumnodes);
         for(i=0;i<vnumnodes-1;i++){
-                node_list[vnumnodes+pnumnodes+i] = pentabase->GetNode(i);
-                cs_list[vnumnodes+pnumnodes+i]   = XYEnum;
+                node_list[i] = pentabase->GetNode(i);
+                cs_list[i]   = XYEnum;
         }
         for(i=0;i<vnumnodes;i++){
-                node_list[i] = element->GetNode(i);
-                cs_list[i]   = XYZEnum;
+                node_list[i+vnumnodes-1] = element->GetNode(i);
+                cs_list[i+vnumnodes-1]   = XYZEnum;
         }
         for(i=0;i<pnumnodes;i++){
-                node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
-                cs_list[vnumnodes+i]   = PressureEnum;
+                node_list[2*vnumnodes-1+i] = element->GetNode(vnumnodes+i);
+                cs_list[2*vnumnodes-1+i]   = PressureEnum;
         }
 
         /*Initialize Element matrix and return if necessary*/
-        ElementMatrix* Ke1=element->NewElementMatrix(FSvelocityEnum);
-        ElementMatrix* Ke2=pentabase->NewElementMatrixCoupling(6,SSAApproximationEnum);
+        ElementMatrix* Ke1=pentabase->NewElementMatrixCoupling(6,SSAApproximationEnum);
+        ElementMatrix* Ke2=element->NewElementMatrix(FSvelocityEnum);
         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
         delete Ke1; delete Ke2;
 
@@ -3698 +3704 @@
                                         &Ke_gg2[0][0],1);
 
         } 
-        for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[i*numdoftotal+j+numdofs]+=Ke_gg[i][j];
-        for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[(i+numdofs)*numdoftotal+j]+=Ke_gg2[i][j];
+        for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
+        for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j];
 
         /*Transform Coordinate System*/
         element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);