Changeset 20658

Timestamp:

05/27/16 14:40:18 (9 years ago)

Author:

Mathieu Morlighem

Message:

CHG: much more general PenaltyCreateKMatrix

File:

• issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp (modified) (2 diffs)

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

@@ -385 +385 @@
 ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceSSAHO(IssmDouble kmax){/*{{{*/
 
-        const int numdof=NUMVERTICES*NDOF2;
+        int        numdof,numdof2,N;
         IssmDouble penalty_offset;
 
@@ -394 +394 @@
         parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum);
 
-        //Create elementary matrix: add penalty to 
-        Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset);
-        Ke->values[0*numdof+2]=-kmax*pow(10.,penalty_offset);
-        Ke->values[2*numdof+0]=-kmax*pow(10.,penalty_offset);
-        Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset);
-
-        Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset);
-        Ke->values[1*numdof+3]=-kmax*pow(10.,penalty_offset);
-        Ke->values[3*numdof+1]=-kmax*pow(10.,penalty_offset);
-        Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset);
+        /*Get number of dof for these two nodes*/
+        numdof =this->nodes[0]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+        numdof2=this->nodes[1]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+        N=NUMVERTICES*numdof;
+
+        /*Add penalty to Element matrix*/
+        for(int i=0;i<numdof;i++){
+                Ke->values[         i*N+i       ]=+kmax*pow(10.,penalty_offset);
+                Ke->values[         i*N+numdof+i]=-kmax*pow(10.,penalty_offset);
+                Ke->values[(numdof+i)*N+i       ]=-kmax*pow(10.,penalty_offset);
+                Ke->values[(numdof+i)*N+numdof+i]=+kmax*pow(10.,penalty_offset);
+        }
 
         /*Clean up and return*/