Changeset 25936

Timestamp:

01/13/21 15:09:59 (4 years ago)

Author:

Cheng Gong

Message:

BUG: fixed a bug in the Nitsche method, where the viscosity should not be multiplied by 2 in the assembly. Also, some improvment in the coding of GLS, but nothing really changed for GLS.

File:

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

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

@@ -4140 +4140 @@
                 mk = 1.0 /3.0;
                 element->ElementSizes(&hx,&hy,&hz);
-                hk = max(max(hx, hy), hz);
-                // if(!element->IsOnDirichletBoundary()) {
-                        Tau = - mk * hk * hk * 0.125 / viscosity;
-
+                // hk = max(max(hx, hy), hz);
+                hk = max(hx, hy);
+                Tau = - mk * hk * hk * 0.125 / viscosity;
 
                 for (int q=0;q<dim;q++){
@@ -4149 +4148 @@
                         for (int p=0;p<dim;p++){
                                 for(int i=0;i<vnumnodes;i++){
-                                        SU[q*vnumnodes*(dim+1)+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
+                                        SU[q*numdof+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
                                 }
                         }
@@ -4155 +4154 @@
                         for(int i=0;i<vnumnodes;i++){
                                 for (int p=0;p<dim;p++){
-                                        SU[q*vnumnodes*(dim+1)+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
+                                        SU[q*numdof+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
                                 }
                         }
                         /* column 4 for the pressure */
-                        for(int i=0;i<vnumnodes;i++){
-                                        SU[q*vnumnodes*(dim+1)+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
+                        for(int i=0;i<pnumnodes;i++){
+                                        SU[q*numdof+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
                         }
                 }
@@ -4226 +4225 @@
         int vnumnodes = element->NumberofNodesVelocity();
         int pnumnodes = element->NumberofNodesPressure();
+        int numdof    = vnumnodes*dim + pnumnodes;
 
         /*Prepare coordinate system list*/
@@ -4291 +4291 @@
                 mk = 1.0 /3.0;
                 element->ElementSizes(&hx,&hy,&hz);
-                hk = max(max(hx, hy), hz);
-                // if(!element->IsOnDirichletBoundary()) {
-                        Tau = - mk * hk * hk * 0.125 / viscosity;
-                // }
-                // else {
-                //      Tau = 0.0;
-                // }
+                // hk = max(max(hx, hy), hz);
+                hk = max(hx, hy);
+                Tau = - mk * hk * hk * 0.125 / viscosity;
 
                 for (int q=0;q<dim;q++){
@@ -4303 +4299 @@
                         for (int p=0;p<dim;p++){
                                 for(int i=0;i<vnumnodes;i++){
-                                        SU[q*vnumnodes*(dim+1)+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
+                                        SU[q*numdof+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
                                 }
                         }
@@ -4309 +4305 @@
                         for(int i=0;i<vnumnodes;i++){
                                 for (int p=0;p<dim;p++){
-                                        SU[q*vnumnodes*(dim+1)+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
+                                        SU[q*numdof+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
                                 }
                         }
                         /* column 4 for the pressure */
-                        for(int i=0;i<vnumnodes;i++){
-                                        SU[q*vnumnodes*(dim+1)+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
+                        for(int i=0;i<pnumnodes;i++){
+                                        SU[q*numdof+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
                         }
                 }
-
                 for(i=0;i<vnumnodes;i++){
                         pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
@@ -4326 +4321 @@
                 }
 
-                for(int i=0;i<vnumnodes*(dim+1);i++){
-                        pe->values[i] += Tau*Jdet*gauss->weight*rho_ice*(-gravity)*SU[(dim-1)*(dim+1)*vnumnodes+i];
+                for(int i=0;i<numdof;i++){
+                        pe->values[i] += Tau*Jdet*gauss->weight*rho_ice*(-gravity)*SU[(dim-1)*numdof+i];
                 }
         }
@@ -5159 +5154 @@
                 for(int i=0;i<vnumnodes;i++){
                         for(int d=0;d<dim;d++){
-                                nsigma[d*vnumnodes+i] = 2.0 * viscosity * bed_normal[d]*vdbasisn[i];
+                                nsigma[d*vnumnodes+i] = 1.0 * viscosity * bed_normal[d]*vdbasisn[i];
                         }
                 }
@@ -5172 +5167 @@
                                         for (int q=0; q<dim; q++) {
                                                 /* gamma*(n*u)*(n*v) */
-                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (-gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i];
+                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i];
                                                 /* sigma(u)*(n*v) */
                                                 Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* nsigma[p*vnumnodes+i] * bed_normal[q] * vbasis[j];