source: issm/oecreview/Archive/15392-16133/ISSM-15474-15475.diff

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

@@ -6917 +6917 @@
         GaussPenta *gauss=NULL;
 
         /*Stabilization*/
-        bool       stabilization = false;
+        bool       stabilization = true;
         IssmDouble dbasis[3][6];
         IssmDouble dmu[3];
         IssmDouble mu;
@@ -6943 +6943 @@
         rigidity=material->GetB();
         diameter=MinEdgeLength(xyz_list);
 
+        if(stabilization){
+                gauss=new GaussPenta();
+                for(int iv=0;iv<6;iv++){
+                        gauss->GaussVertex(iv);
+                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
+                        for(i=0;i<6;i++){
+                                for(j=0;j<3;j++){
+                                        dnodalbasis[i][iv][j] = dbasis[j][i];
+                                }
+                        }
+                }
+                delete gauss;
+        }
+
         /* Start  looping on the number of gaussian points: */
         gauss=new GaussPenta(3,2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -6971 +6985 @@
                 if(stabilization){
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
                         dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
-                        mu = viscosity;
-                        for(i=0;i<6;i++){
-                                for(p=0;p<6;p++){
-                                        for(j=0;j<3;j++){
-                                                dnodalbasis[i][p][j] = dbasis[j][i];
-                                        }
-                                }
-                        }
+                        mu = 2.*viscosity;
                         for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
                                 SU[p][i][j]=0.;
                                 SW[p][i][j]=0.;
@@ -7017 +7024 @@
                         }
 
                 }
-
-
-                //viscosity = 1000*rigidity;
-                //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
-                //if(id==24){
-                //printf("tau %g\n",1./3.*pow(diameter,2)/(8.*2.*viscosity));
-                //}
-                //GetBConduct(&B_stab[0][0],&xyz_list[0][0],gauss); 
-
-                //D_stab[0][0]=D_scalar_stab; D_stab[0][1]=0;             D_stab[0][2]=0;
-                //D_stab[1][0]=0;             D_stab[1][1]=D_scalar_stab; D_stab[1][2]=0;
-                //D_stab[2][0]=0;             D_stab[2][1]=0;             D_stab[2][2]=D_scalar_stab;
-
-                //TripleMultiply(&B_stab[0][0],3,6,1,
-                //                      &D_stab[0][0],3,3,0,
-                //                      &B_stab[0][0],3,6,0,
-                //                      &Ke_temp_stab[0][0],0);
-
-                //for(i=0;i<NUMVERTICES;i++) for(j=0;j<NUMVERTICES;j++) Ke->values[numdof*(i*4+3)+j*4+3]+=Ke_temp_stab[i][j];
-
         }
 
-//      if(id==24){
-//              _error_("");
-//      }
-
-        //Ke->Echo();
-        //_error_("stop");
         /*Transform Coordinate System*/
         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
         if(id==24){
@@ -7943 +7924 @@
         /*Find minimal length*/
         diameter=MinEdgeLength(xyz_list);
 
+        if(stabilization){
+                gauss=new GaussPenta();
+                for(int iv=0;iv<6;iv++){
+                        gauss->GaussVertex(iv);
+                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
+                        for(i=0;i<6;i++){
+                                for(j=0;j<3;j++){
+                                        dnodalbasis[i][iv][j] = dbasis[j][i];
+                                }
+                        }
+                }
+                delete gauss;
+        }
+
         /* Start  looping on the number of gaussian points: */
-        gauss=new GaussPenta(2,3);
+        gauss=new GaussPenta(3,2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 
                 gauss->GaussPoint(ig);
@@ -7969 +7964 @@
                 if(stabilization){
                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
                         dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
-                        mu = viscosity;
-                        for(i=0;i<6;i++){
-                                for(p=0;p<6;p++){
-                                        for(j=0;j<3;j++){
-                                                dnodalbasis[i][p][j] = dbasis[j][i];
-                                        }
-                                }
-                        }
-                        //dNodalBasisdx(1:n,p,:) = dBasisdx(1:n,:)
+                        mu = 2.*viscosity;
                         for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
                                 SW[p][i][j]=0.;
                         }
@@ -7988 +7975 @@
                                                 SW[p][i][i] += -dmu[j]*dbasis[j][p];
                                                 SW[p][j][i] += -dmu[j]*dbasis[i][p];
                                                 for(ii=0;ii<6;ii++){
-                                                        //SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
-                                                        //SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
-                                                        SW[p][i][i] += -mu*dnodalbasis[p][ii][j];
-                                                        SW[p][j][i] += -mu*dnodalbasis[p][ii][i];
+                                                        SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
+                                                        SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
                                                 }
                                         }
                                 }
@@ -8006 +7991 @@
                         }
 
                 }
-                ///*Add stabilization*/
-                //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
-                //GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0],gauss);
-
-                //for(i=0;i<NUMVERTICES;i++){
-                //      Pe_gaussian[i*NDOF4+3]+=-rho_ice*gravity*D_scalar_stab*dh1dh6[2][i];
-                //      Pe_gaussian[i*NDOF4+3]+=forcex*D_scalar_stab*dh1dh6[0][i];
-                //      Pe_gaussian[i*NDOF4+3]+=forcey*D_scalar_stab*dh1dh6[1][i];
-                //      Pe_gaussian[i*NDOF4+3]+=forcez*D_scalar_stab*dh1dh6[2][i];
-                //}
         }
 
         /*Transform coordinate system*/
         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
 
-        if(id==24){
-                printarray(pe->values,24,1);
-                _error_("STOP");
-        }
-
         /*Clean up and return*/
         delete gauss;
         return pe;