Changeset 24120

Timestamp:

08/23/19 17:04:24 (6 years ago)

Author:

tsantos

Message:

BUG: fixed matrix computation for Streamline Upwind stabilization scheme

File:

• issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (modified) (3 diffs)

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

@@ -376 +376 @@
                                 break;
                         case 2:
+                                /*Streamline upwinding*/
+                                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+                                vxaverage_input->GetInputAverage(&vx);
                                 if(dim==1){
                                         vel=fabs(vx)+1.e-8;
-                                        D[0]=h/(2*vel)*vx*vx;
                                 }
                                 else{
-                                        /*Streamline upwinding*/
+                                        vyaverage_input->GetInputAverage(&vy);
                                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
-                                        D[0*dim+0]=h/(2*vel)*vx*vx;
-                                        D[1*dim+0]=h/(2*vel)*vy*vx;
-                                        D[0*dim+1]=h/(2*vel)*vx*vy;
-                                        D[1*dim+1]=h/(2*vel)*vy*vy;
                                 }
+                                tau=h/(2*vel);
                                 break;
                         case 5:
@@ -403 +402 @@
                                 _error_("Stabilization "<<stabilization<<" not supported yet");
                 }
-                if(stabilization==1 || stabilization==2){
+                if(stabilization==1){
+                        /*SSA*/
                         if(dim==1) D[0]=D_scalar*D[0];
                         else{
@@ -416 +416 @@
                                                 Bprime,dim,numnodes,0,
                                                 &Ke->values[0],1);
+                }
+                if(stabilization==2){
+                        /*Streamline upwind*/
+                        for(int i=0;i<numnodes;i++){
+                                for(int j=0;j<numnodes;j++){
+                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
+                                }
+                        }
                 }
                 if(stabilization==5){/*{{{*/