Changeset 602

Timestamp:

05/26/09 17:12:55 (15 years ago)

Author:

seroussi

Message:

fixed a few problems with prognostic

Location:

issm/trunk/src/c/objects

Files:

• Tria.cpp (modified) (8 diffs)
• Tria.h (modified) (1 diff)

issm/trunk/src/c/objects/Tria.cpp

@@ -551 +551 @@
         /* matrices: */
         double L[numgrids];
+        double B[2][numgrids];
+        double Bprime[2][numgrids];
         double DL[2][2]={0.0};
         double DLprime[2][2]={0.0};
@@ -556 +558 @@
         double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix 
         double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
-        double Ke_gg_thickness[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+        double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+        double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
         
         double Jdettria;
         
         /*input parameters for structural analysis (diagnostic): */
-        double  vxvy_list[numgrids][2]={0,0};
-        double  vx_list[numgrids]={0,0};
-        double  vy_list[numgrids]={0,0};
+        double  vxvy_list[numgrids][2]={0.0};
+        double  vx_list[numgrids]={0.0};
+        double  vy_list[numgrids]={0.0};
         double  dvx[2];
         double  dvy[2];
         double  vx,vy;
+        double  dvxdx,dvydy;
         double  v_gauss[2]={0.0};
         double  K[2][2]={0.0};
@@ -636 +640 @@
                 
                 /*Get B  and B prime matrix: */
-                //GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
-                //GetBprime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
+                GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
+                GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
                 
                 //Get vx, vy and their derivatives at gauss point
@@ -646 +650 @@
                 GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
 
-                //dvxdx=dvx[0];
-                //dvydy=dvy[1];
+                dvxdx=dvx[0];
+                dvydy=dvy[1];
 
                 DL_scalar=dt*gauss_weight*Jdettria;
 
                 //Create DL and DLprime matrix
-                //DL[0][0]=DL_scalar*dvxdx;
-                //DL[1][1]=DL_scalar*dvydy;
+                DL[0][0]=DL_scalar*dvxdx;
+                DL[1][1]=DL_scalar*dvydy;
 
                 DLprime[0][0]=DL_scalar*vx;
@@ -661 +665 @@
                 //Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;
 
-                /*TripleMultiply( &B[0][0],numdof,numdof,1,
-                                          &DL_scalar,1,1,0,
-                                          &L[0],1,numdof,0,
-                                          &Ke_gg_gaussian[0][0],0);*/
-
-
-
-
+                TripleMultiply( &B[0][0],2,numdof,1,
+                                          &DL[0][0],2,2,0,
+                                          &B[0][0],2,numdof,0,
+                                          &Ke_gg_thickness1[0][0],0);
+
+                TripleMultiply( &B[0][0],2,numdof,1,
+                                          &DLprime[0][0],2,2,0,
+                                          &Bprime[0][0],2,numdof,0,
+                                          &Ke_gg_thickness2[0][0],0);
 
                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
+                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
                 
                 #ifdef _DEBUGELEMENTS_
@@ -696 +703 @@
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-
-        /*Do not forget to include friction: */
-        if(!shelf){
-                //CreateKMatrixPrognosticFriction(Kgg,inputs,analysis_type,sub_analysis_type);
-        }
-
-
         #ifdef _DEBUGELEMENTS_
         if(my_rank==RANK && id==ELID){ 
@@ -766 +766 @@
         double  thickness_g;
         double  dt;
-        int     dofs[1]={1};
+        int     dofs[1]={0};
         int     found=0;
 
@@ -1706 +1706 @@
         }
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetB_prog"
+
+void Tria::GetB_prog(double* B_prog, double* xyz_list, double* gauss_l1l2l3){
+
+        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+         * For grid i, Bi can be expressed in the basic coordinate system
+         * by: 
+         *       Bi_basic=[ h ]
+         *                [ h ]
+         * where h is the interpolation function for grid i.
+         *
+         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
+         */
+        
+        int i;
+        const int NDOF1=1;
+        const int numgrids=3;
+
+        double l1l2l3[numgrids];
+
+
+        /*Get dh1dh2dh3 in basic coordinate system: */
+        GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
+
+        #ifdef _DEBUG_ 
+        for (i=0;i<3;i++){
+                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
+        }
+        #endif
+
+        /*Build B_prog: */
+        for (i=0;i<numgrids;i++){
+                *(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
+                *(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
+        }
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetBPrime_prog"
+
+void Tria::GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3){
+
+        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+         * For grid i, Bi' can be expressed in the basic coordinate system
+         * by: 
+         *       Bi_prime__basic=[ dh/dx ]
+         *                       [ dh/dy ]
+         * where h is the interpolation function for grid i.
+         *
+         * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
+         */
+        
+        int i;
+        const int NDOF1=1;
+        const int NDOF2=2;
+        const int numgrids=3;
+
+        /*Same thing in the basic coordinate system: */
+        double dh1dh2dh3_basic[NDOF2][numgrids];
+
+        /*Get dh1dh2dh3 in basic coordinates system : */
+        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list,gauss_l1l2l3);
+
+        /*Build B': */
+        for (i=0;i<numgrids;i++){
+                *(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh2dh3_basic[0][i]; 
+                *(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh2dh3_basic[1][i]; 
+        }
+}
+
 
 #undef __FUNCT__ 

issm/trunk/src/c/objects/Tria.h

@@ -84 +84 @@
                 void  GetBPrime(double* Bprime, double* xyz_list, double* gauss_l1l2l3);
                 void  GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof);
+                void  GetB_prog(double* B_prog, double* xyz_list, double* gauss_l1l2l3);
+                void  GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3);
                 void  GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3);
                 void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3);