Changeset 10527

Timestamp:

11/08/11 11:46:01 (13 years ago)

Author:

Mathieu Morlighem

Message:

Extended CS rotation to macayeal Stokes coupling

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (10 diffs)

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -5346 +5346 @@
 
         /*Constants*/
-        const int    numdofm=NDOF2*NUMVERTICES2D;
-        const int    numdofs=NDOF4*NUMVERTICES;
-        const int    numdoftotal=2*numdofm+numdofs;
+        const int numnodes    = 2 *NUMVERTICES;
+        const int numdofm     = NDOF2 *NUMVERTICES2D;
+        const int numdofs     = NDOF4 *NUMVERTICES;
+        const int numdoftotal = 2 *numdofm+numdofs;
 
         /*Intermediaries */
@@ -5369 +5370 @@
         GaussPenta *gauss=NULL;
         GaussTria  *gauss_tria=NULL;
+        Node       *node_list[numnodes];
+        int         cs_list[numnodes];
 
         /*Find penta on bed as stokes must be coupled to the dofs on the bed: */
         Penta* pentabase=GetBasalElement();
         Tria* tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
+
+        /*Prepare node list*/
+        for(i=0;i<NUMVERTICES;i++){
+                node_list[i+0*NUMVERTICES] = pentabase->nodes[i];
+                node_list[i+1*NUMVERTICES] = this->nodes[i];
+                cs_list[i+0*NUMVERTICES] = XYEnum;
+                cs_list[i+1*NUMVERTICES] = XYZPEnum;
+        }
 
         /*Initialize Element matrix and return if necessary*/
@@ -5425 +5436 @@
         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*/
+        TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+
         /*Clean-up and return*/
         delete tria->matice; delete tria;
@@ -5436 +5450 @@
 
         /*Constants*/
-        const int numdof=NUMVERTICES*NDOF4;
-        const int numdofm=NUMVERTICES*NDOF2;
-        const int numdof2d=NUMVERTICES2D*NDOF4;
-        const int numdof2dm=NUMVERTICES2D*NDOF2;
-        const int numdoftot=numdof+numdofm;
+        const int numnodes  = 2 *NUMVERTICES;
+        const int numdof    = NUMVERTICES *NDOF4;
+        const int numdofm   = NUMVERTICES *NDOF2;
+        const int numdof2d  = NUMVERTICES2D *NDOF4;
+        const int numdof2dm = NUMVERTICES2D *NDOF2;
+        const int numdoftot = numdof+numdofm;
 
         /*Intermediaries */
@@ -5461 +5476 @@
         Friction*  friction=NULL;
         GaussPenta *gauss=NULL;
+        Node       *node_list[numnodes];
+        int         cs_list[numnodes];
 
         /*If on water or not Stokes, skip stiffness: */
@@ -5469 +5486 @@
         ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
         delete Ke1; delete Ke2;
+
+        /*Prepare node list*/
+        for(i=0;i<NUMVERTICES;i++){
+                node_list[i+0*NUMVERTICES] = this->nodes[i];
+                node_list[i+1*NUMVERTICES] = this->nodes[i];
+                cs_list[i+0*NUMVERTICES] = XYEnum;
+                cs_list[i+1*NUMVERTICES] = XYZPEnum;
+        }
 
         /*Retrieve all inputs and parameters*/
@@ -5527 +5552 @@
                 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2dm;j++) Ke->values[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j];
         }
+
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
 
         /*Clean up and return*/
@@ -6324 +6352 @@
         }
 
+        /*Transform coordinate system*/
+        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
+
         /*Clean up and return*/
         delete gauss;
@@ -6394 +6425 @@
                 }
         }
+
+        /*Transform coordinate system*/
+        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
 
         /*Clean up and return*/
@@ -7474 +7508 @@
         }
 
+        /*Transform solution in Cartesian Space*/
+        TransformSolutionCoord(&macayeal_values[0],this->nodes,NUMVERTICES,XYEnum);
+        TransformSolutionCoord(&stokes_values[0],this->nodes,NUMVERTICES,XYZPEnum);
+
         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
         for(i=0;i<NUMVERTICES;i++){