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

Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16908)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp       (revision 16909)
@@ -3089,8 +3089,82 @@
        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,7 +3183,6 @@
        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,9 +3233,7 @@
                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,7 +3247,7 @@
        delete gauss;
        delete gauss_tria;
        xDelete<IssmDouble>(Ke_gg);
-       xDelete<IssmDouble>(Ke_gg_gaussian);
+       xDelete<IssmDouble>(xyz_list);
        xDelete<int>(cs_list);
        return Ke;