Changeset 10524

Timestamp:

11/08/11 11:23:34 (13 years ago)

Author:

Mathieu Morlighem

Message:

Extended CS rotation to macayeal pattun coupling

Location:

issm/trunk/src/c

Files:

• objects/Elements/Penta.cpp (modified) (8 diffs)
• shared/Elements/CoordinateSystemTransform.cpp (modified) (1 diff)

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

@@ -5141 +5141 @@
 
         /*Constants*/
-        const int    numdofm=NDOF2*NUMVERTICES2D;
-        const int    numdofp=NDOF2*NUMVERTICES;
-        const int    numdoftotal=2*NDOF2*NUMVERTICES;
+        const int numnodes    = 2 *NUMVERTICES;
+        const int numdofm     = NDOF2 *NUMVERTICES2D;
+        const int numdofp     = NDOF2 *NUMVERTICES;
+        const int numdoftotal = 2 *NDOF2*NUMVERTICES;
 
         /*Intermediaries */
@@ -5159 +5160 @@
         GaussPenta *gauss=NULL;
         GaussTria  *gauss_tria=NULL;
+        Node       *node_list[numnodes];
+        int         cs_list[numnodes];
 
         /*Find penta on bed as pattyn 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.
+        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] = XYEnum;
+        }
 
         /*Initialize Element matrix*/
@@ -5209 +5220 @@
         for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j][i];
 
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+
         /*Clean-up and return*/
         delete tria->matice; delete tria;
@@ -5220 +5234 @@
 
         /*Constants*/
-        const int numdof        = NDOF2 *NUMVERTICES;
-        const int numdoftotal   = NDOF4 *NUMVERTICES;
+        const int numnodes    = 2 *NUMVERTICES;
+        const int numdof      = NDOF2 *NUMVERTICES;
+        const int numdoftotal = NDOF4 *NUMVERTICES;
         
         /*Intermediaries */
@@ -5238 +5253 @@
         Friction  *friction = NULL;
         GaussPenta *gauss=NULL;
+        Node       *node_list[numnodes];
+        int         cs_list[numnodes];
 
         /*Initialize Element matrix and return if necessary*/
@@ -5245 +5262 @@
         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] = XYEnum;
+        }
 
         /*retrieve inputs :*/
@@ -5293 +5318 @@
         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j];
         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j];
+
+        /*Transform Coordinate System*/
+        TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
 
         /*Clean up and return*/
@@ -7354 +7382 @@
         }
 
+        /*Transform solution in Cartesian Space*/
+        TransformSolutionCoord(&macayeal_values[0],penta->nodes,NUMVERTICES,XYEnum);
+        TransformSolutionCoord(&pattyn_values[0],   this->nodes,NUMVERTICES,XYEnum);
+
         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
         for(i=0;i<NUMVERTICES;i++){

issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp

@@ -45 +45 @@
                         case XYEnum:
                                 /*We remove the z component, we need to renormalize x and y*/
-                                x_norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]);
-                                y_norm = sqrt( coord_system[0][1]*coord_system[0][1] + coord_system[1][1]*coord_system[1][1]);
+                                x_norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]); _assert_(x_norm>1.e-5);
+                                y_norm = sqrt( coord_system[0][1]*coord_system[0][1] + coord_system[1][1]*coord_system[1][1]); _assert_(y_norm>1.e-5);
                                 transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0]/x_norm;
                                 transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1]/y_norm;