Changeset 18241

Timestamp:

07/14/14 14:07:50 (11 years ago)

Author:

seroussi

Message:

NEW: preparing friction for LA

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (7 diffs)

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

@@ -7 +7 @@
 #include "../cores/cores.h"
 
-//#define FSANALYTICAL 24
+//#define FSANALYTICAL 12
 
 /*Model processing*/
@@ -3246 +3246 @@
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
-        IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
-        IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+        IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
+        IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
 
         /*Retrieve all inputs and parameters*/
@@ -3271 +3271 @@
                 this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-                for(int i=0;i<dim-1;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
-
-                TripleMultiply(B,dim-1,numdof,1,
-                                        D,dim-1,dim-1,0,
-                                        B,dim-1,numdof,0,
+                for(int i=0;i<dim;i++) D[i*dim+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
+
+                TripleMultiply(B,dim,numdof,1,
+                                        D,dim,dim,0,
+                                        B,dim,numdof,0,
                                         &Ke->values[0],1);
         }
@@ -3394 +3394 @@
         /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
-        IssmDouble*    B  = xNew<IssmDouble>((dim-1)*numdof);
-        IssmDouble*    D  = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+        IssmDouble*    B  = xNew<IssmDouble>(dim*numdof);
+        IssmDouble*    D  = xNewZeroInit<IssmDouble>(dim*dim);
 
         /*Retrieve all inputs and parameters*/
@@ -3431 +3431 @@
                 this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-                for(int i=0;i<dim-1;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
-
-                TripleMultiply(B,dim-1,numdof,1,
-                                        D,dim-1,dim-1,0,
-                                        B,dim-1,numdof,0,
+                for(int i=0;i<dim;i++) D[i*dim+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
+
+                TripleMultiply(B,dim,numdof,1,
+                                        D,dim,dim,0,
+                                        B,dim,numdof,0,
                                         &Ke->values[0],1);
         }
@@ -4429 +4429 @@
                         B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
                         B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
+
+                        B[(vnumdof+pnumdof)*2+3*i+0] = 0.;
+                        B[(vnumdof+pnumdof)*2+3*i+1] = 0.;
+                        B[(vnumdof+pnumdof)*2+3*i+2] = vbasis[i];
                 }
                 for(int i=0;i<pnumnodes;i++){
@@ -4437 +4441 @@
         else{
                 for(int i=0;i<vnumnodes;i++){
-                        B[2*i+0] = vbasis[i];
-                        B[2*i+1] = 0.;
+                        B[(vnumdof+pnumdof)*0+2*i+0] = vbasis[i];
+                        B[(vnumdof+pnumdof)*0+2*i+1] = 0.;
+
+                        B[(vnumdof+pnumdof)*1+2*i+0] = 0.;
+                        B[(vnumdof+pnumdof)*1+2*i+1] = vbasis[i];
                 }