Changeset 16373

Timestamp:

10/10/13 16:39:59 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added friction for FS flowline

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Penta.cpp (modified) (3 diffs)
• Elements/Tria.cpp (modified) (4 diffs)
• Elements/TriaRef.cpp (modified) (1 diff)
• Elements/TriaRef.h (modified) (1 diff)
• Loads/Friction.cpp (modified) (1 diff)

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

@@ -8186 +8186 @@
         int         analysis_type,approximation;
         IssmDouble  alpha2,Jdet2d;
-        IssmDouble  FSreconditioning,viscosity;
-        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
@@ -8213 +8211 @@
         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -8230 +8227 @@
                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
                 GetLFS(BFriction,gauss);
-
-                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3208 +3208 @@
         /*Intermediaries */
         int         i,j;
-        int         analysis_type,approximation;
+        int         analysis_type;
         int         indices[2];
         IssmDouble  alpha2,Jdet;
-        IssmDouble  FSreconditioning,viscosity;
-        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_seg[NUMVERTICES1D][3];
@@ -3232 +3230 @@
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke        = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
-        IssmDouble*    BFriction = xNew<IssmDouble>(2*numdof);
-        IssmDouble*    D         = xNewZeroInit<IssmDouble>(2*2);
+        IssmDouble*    BFriction = xNew<IssmDouble>(1*numdof);
+        IssmDouble     D;
 
         /*Retrieve all inputs and parameters*/
         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -3256 +3253 @@
 
                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
-
-                _error_("STOP");
+                GetLFS(BFriction,gauss);
+
+                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
+                D = +alpha2*gauss->weight*Jdet; //taub_x = -alpha2 vx
+
+                TripleMultiply(BFriction,1,numdof,1,
+                                        &D,1,1,0,
+                                        BFriction,1,numdof,0,
+                                        Ke->values,1);
         }
 
@@ -3267 +3271 @@
         delete friction;
         xDelete<IssmDouble>(BFriction);
-        xDelete<IssmDouble>(D);
         xDelete<int>(cs_list);
         return Ke;

issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp

@@ -469 +469 @@
         /*Clean-up*/
         xDelete<IssmDouble>(basis);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetLFS{{{*/
+void TriaRef::GetLFS(IssmDouble* LFS, GaussTria* gauss){
+        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+         * For node i, Li can be expressed in the actual coordinate system
+         * by: 
+         *       Li=[ h 0 0 ]
+         * where h is the interpolation function for node i.
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int pnumnodes = this->NumberofNodesPressure();
+        int vnumnodes = this->NumberofNodesVelocity();
+        int pnumdof   = pnumnodes;
+        int vnumdof   = vnumnodes*NDOF2;
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
+        GetNodalFunctionsVelocity(vbasis,gauss);
+
+        /*Build LFS: */
+        for(int i=0;i<vnumnodes;i++){
+                LFS[2*i+0] = vbasis[i];
+                LFS[2*i+1] = 0.;
+        }
+
+        for(int i=0;i<pnumnodes;i++){
+                LFS[i+vnumdof+0] = 0.;
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(vbasis);
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/TriaRef.h

@@ -37 +37 @@
                 void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss);
                 void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
+                void GetLFS(IssmDouble* LFS, GaussTria* gauss);
                 void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
                 void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss,int finiteelement);

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -91 +91 @@
         if (Neff<0)Neff=0;
 
-        if(strcmp(element_type,"2d")==0){
+        if(strcmp(element_type,"1d")==0){
+                this->GetInputValue(&vx, gauss,vxenum);
+                vmag=sqrt(vx*vx);
+        }
+        else if(strcmp(element_type,"2d")==0){
                 this->GetInputValue(&vx, gauss,vxenum);
                 this->GetInputValue(&vy, gauss,vyenum);