Changeset 26141

Timestamp:

03/23/21 15:24:26 (4 years ago)

Author:

tsantos

Message:

CHG: working on mono layer HO

File:

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

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

@@ -2735 +2735 @@
         for(int i=0;i<numnodes;i++){
       for(int j=0;j<numnodes;j++){
-         Ke->values[(i+0)*2*2*numnodes+j+0] = KeSSA->values[2*i*2*numnodes+2*j]; //K11
-         Ke->values[(i+6)*2*2*numnodes+j+6] = KeSSA->values[(2*i+1)*2*numnodes+2*j+1]; //K33
+         Ke->values[(i+0*numnodes)*2*2*numnodes+j+0*numnodes] = KeSSA->values[2*i*2*numnodes+2*j]; //K11
+         Ke->values[(i+2*numnodes)*2*2*numnodes+j+2*numnodes] = KeSSA->values[(2*i+1)*2*numnodes+2*j+1]; //K33
       }  
    }
@@ -2787 +2787 @@
                 for(int i=0;i<numnodes;i++){//shape functions on tria element
                         for(int j=0;j<numnodes;j++){
-                                Ke->values[(i+0)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+0*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         );//K11
-                                Ke->values[(i+0)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+0*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*(n+1)/(n+2);//K12
-                                Ke->values[(i+0)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+0*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                                                         );//K13
-                                Ke->values[(i+0)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+0*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                      2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                      )*(n+1)/(n+2);//K14
                                 
-                                Ke->values[(i+3)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+1*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*(n+1)/(n+2);//K21
-                                Ke->values[(i+3)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+1*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*2*pow(n+1,2)/((2*n+3)*(n+2)) 
@@ -2809 +2809 @@
                                                         gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
                                                         ;//K22
-                                Ke->values[(i+3)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+1*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*(n+1)/(n+2);//K23
-                                Ke->values[(i+3)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+1*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
                                                         )*2*pow(n+1,2)/((2*n+3)*(n+2));//K24
                                 
-                                Ke->values[(i+6)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+2*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                                                         );//K31
-                                Ke->values[(i+6)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+2*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*(n+1)/(n+2);//K32
-                                Ke->values[(i+6)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+2*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         );//K33
-                                Ke->values[(i+6)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+2*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*(n+1)/(n+2);//K34
 
-                                Ke->values[(i+9)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+3*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                      2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                      )*(n+1)/(n+2);//K41
-                                Ke->values[(i+9)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+3*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                      2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
                      )*2*pow(n+1,2)/((2*n+3)*(n+2));//K42
-                                Ke->values[(i+9)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+3*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*(n+1)/(n+2);//K43
-                                Ke->values[(i+9)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
+                                Ke->values[(i+3*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(
                                                         4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
                                                         )*2*pow(n+1,2)/((2*n+3)*(n+2))
@@ -2926 +2926 @@
 
                 for(int i=0;i<numnodes;i++){
-                        pe->values[i+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1
-                        pe->values[i+3]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F2
-                        pe->values[i+6]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; //F3
-                        pe->values[i+9]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F4
+                        pe->values[i+0*numnodes]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1
+                        pe->values[i+1*numnodes]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F2
+                        pe->values[i+2*numnodes]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; //F3
+                        pe->values[i+3*numnodes]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F4
                 }
         }
@@ -2999 +2999 @@
 
                 for (int i=0;i<numnodes;i++){
-                        pe->values[i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; // F1
-                        pe->values[i+3]+= pressure_sh*Jdet*gauss->weight*normal[0]*basis[i]; // F2
-                        pe->values[i+6]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; // F3
-                        pe->values[i+9]+= pressure_sh*Jdet*gauss->weight*normal[1]*basis[i]; // F4
+                        pe->values[i+0*numnodes]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; // F1
+                        pe->values[i+1*numnodes]+= pressure_sh*Jdet*gauss->weight*normal[0]*basis[i]; // F2
+                        pe->values[i+2*numnodes]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; // F3
+                        pe->values[i+3*numnodes]+= pressure_sh*Jdet*gauss->weight*normal[1]*basis[i]; // F4
                 }
         }