source: issm/oecreview/Archive/16554-17801/ISSM-16908-16909.diff@ 17802

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -3089 +3089 @@
         return Ke;
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/
-        return NULL;
 
+        if(element->IsFloating() || !element->IsOnBed()) return NULL;
+
+        /*Constants*/
+        int numnodes    = element->GetNumberOfNodes();
+        int numdof      = 2*numnodes;
+        int numdoftotal = 4*numnodes;
+
+        /*Intermediaries */
+        int         i,j;
+        IssmDouble  Jdet2d,alpha2;
+        IssmDouble *xyz_list_tria=NULL;
+        IssmDouble* L           = xNewZeroInit<IssmDouble>(2*numdof);
+        IssmDouble  DL[2][2]     ={{ 0,0 },{0,0}}; //for basal drag
+        IssmDouble  DL_scalar;
+        IssmDouble* Ke_gg       = xNewZeroInit<IssmDouble>(numdof*numdof);
+        Node       *node_list[2*numnodes];
+        int*        cs_list= xNew<int>(2*numnodes);
+
+        /*Initialize Element matrix and return if necessary*/
+        ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
+        ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
+        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+        delete Ke1; delete Ke2;
+
+        /*Prepare node list*/
+        for(i=0;i<numnodes;i++){
+                node_list[i+0*numnodes] = element->GetNode(i);
+                node_list[i+1*numnodes] = element->GetNode(i);
+                cs_list[i+0*numnodes] = XYEnum;
+                cs_list[i+1*numnodes] = XYEnum;
+        }
+
+        /*retrieve inputs :*/
+        element->GetVerticesCoordinatesBase(&xyz_list_tria);
+        Input* vx_input=element->GetInput(VxEnum);  _assert_(vx_input);
+        Input* vy_input=element->GetInput(VyEnum);  _assert_(vy_input);
+        Input* vz_input=element->GetInput(VzEnum);  _assert_(vz_input);
+
+        /*build friction object, used later on: */
+        Friction* friction=new Friction(element,2);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGaussBase(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /*Friction: */
+                friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+                element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss);
+                this->GetBHOFriction(L,element,xyz_list_tria,gauss);
+
+                DL_scalar=alpha2*gauss->weight*Jdet2d;
+                for (i=0;i<2;i++) DL[i][i]=DL_scalar; 
+
+                /*  Do the triple producte tL*D*L: */
+                TripleMultiply( L,2,numdof,1,
+                                        &DL[0][0],2,2,0,
+                                        L,2,numdof,0,
+                                        Ke_gg,1);
+        }
+
+        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j];
+        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j];
+
+        /*Transform Coordinate System*/
+        element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
+
+        /*Clean up and return*/
+        delete gauss;
+        delete friction;
+        xDelete<int>(cs_list);
+        xDelete<IssmDouble>(xyz_list_tria);
+        xDelete<IssmDouble>(Ke_gg);
+        return Ke;
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/
 
@@ -3109 +3183 @@
         IssmDouble  D[3][3]={0.0}; // material matrix, simple scalar matrix.
         IssmDouble  D_scalar;
         IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
-        IssmDouble* Ke_gg_gaussian = xNew<IssmDouble>(numdofp*numdofm);
         Node        *node_list[2*numnodes];
         int*         cs_list= xNew<int>(2*numnodes);
 
@@ -3160 +3233 @@
                 TripleMultiply( B,3,numdofp,1,
                                         &D[0][0],3,3,0,
                                         Bprime,3,numdofm,0,
-                                        Ke_gg_gaussian,0);
-
-                for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i*numdofm+j]+=Ke_gg_gaussian[i*numdofm+j];
+                                        Ke_gg,1);
         } 
         for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
         for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
@@ -3176 +3247 @@
         delete gauss;
         delete gauss_tria;
         xDelete<IssmDouble>(Ke_gg);
-        xDelete<IssmDouble>(Ke_gg_gaussian);
+        xDelete<IssmDouble>(xyz_list);
         xDelete<int>(cs_list);
         return Ke;