Changeset 15776

Timestamp:

08/09/13 13:24:34 (12 years ago)

Author:

seroussi

Message:

BUG: fixed coupling HO/FS ice sheet

File:

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

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

@@ -6710 +6710 @@
         Tria*  tria=pentabase->SpawnTria(0); //lower face is 0, upper face is 1.
 
-        /*Prepare node list*/
+        /*prepare node list*/
         for(i=0;i<NUMVERTICES;i++){
                 node_list[i+0*NUMVERTICES] = pentabase->nodes[i];
@@ -7112 +7112 @@
 
         /*Constants*/
-        const int numnodes       = 2 *NUMVERTICES;
+        const int numnodes       = 3 *NUMVERTICES+1;
         const int numdofp        = NDOF2 *NUMVERTICES;
         const int numdofs        = NDOF4 * 6 + NDOF3;
@@ -7119 +7119 @@
         /*Intermediaries*/
         int       i,j,init;
+        Node  *node_list[NUMVERTICES*3+1];
+        int   cs_list[NUMVERTICES*3+1];
+        int   cs_list2[NUMVERTICES*2+1];
 
         /*Some parameters needed*/
         init=this->element_type;
+
+        /*prepare node list*/
+        for(i=0;i<NUMVERTICES+1;i++){
+                node_list[i+NUMVERTICES] = this->nodes[i];
+                cs_list[i+NUMVERTICES]   = XYZEnum;
+                cs_list2[i]              = XYZEnum;
+        }
+        for(i=0;i<NUMVERTICES;i++){
+                node_list[i]                 = this->nodes[i];
+                node_list[i+2*NUMVERTICES+1] = this->nodes[i+NUMVERTICES+1];
+                cs_list[i]                   = XYEnum;
+                cs_list[i+2*NUMVERTICES+1]   = PressureEnum;
+                cs_list2[i+NUMVERTICES+1]    = PressureEnum;
+        }
 
         /*compute all stiffness matrices for this element*/
@@ -7131 +7148 @@
         /*Compute HO Matrix with P1 element type\n");*/
         this->element_type=P1Enum;
-        Ke1=CreateKMatrixStressbalanceHO();
+        Ke1=CreateKMatrixStressbalanceHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
         this->element_type=init;
         /*Compute FS Matrix and condense it \n");*/
-        Ke2=CreateKMatrixStressbalanceFS(); 
+        Ke2=CreateKMatrixStressbalanceFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,2*NUMVERTICES+1,cs_list2);
         int indices[3]={18,19,20};
         Ke2->StaticCondensation(3,&indices[0]);
@@ -7148 +7165 @@
 
         /*Transform Coordinate System*/ //Do not transform, already done in the matrices
-        //TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+        TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
 
         /*clean-up and return*/
@@ -8258 +8275 @@
         ElementVector* pe1=CreatePVectorCouplingHOFSViscous();
         ElementVector* pe2=CreatePVectorCouplingHOFSFriction();
-//      if(id==15){
-//      pe1->Echo();
-//      pe2->Echo();
-//      _error_("");
-//      }
         ElementVector* pe =new ElementVector(pe1,pe2);