Changeset 21797

Timestamp:

07/14/17 10:41:25 (8 years ago)

Author:

seroussi

Message:

BUG: confusion between CG and DG

File:

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

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

@@ -469 +469 @@
 
         /*Intermediaries */
+        int         migration_style,point1;
+        bool        mainlyfloating;
+        IssmDouble  fraction1,fraction2;
         IssmDouble  Jdet,dt;
         IssmDouble  ms,mb,gmb,fmb,thickness;
         IssmDouble  gllevelset,phi=1.;
         IssmDouble* xyz_list = NULL;
+        Gauss*      gauss     = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -490 +494 @@
         Input* thickness_input     = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
 
+        /*Recover portion of element that is grounded*/
+        phi=element->GetGroundedPortion(xyz_list);
+        if(migration_style==SubelementMigration2Enum){
+                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+           gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
+        }
+        else{
+                gauss = element->NewGauss(3);
+        }
+
         /* Start  looping on the number of gaussian points: */
-        Gauss* gauss=element->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
@@ -501 +515 @@
                 gmb_input->GetInputValue(&gmb,gauss);
                 fmb_input->GetInputValue(&fmb,gauss);
-                gllevelset_input->GetInputValue(&phi,gauss);
+                gllevelset_input->GetInputValue(&gllevelset,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
-                if(phi>0.) mb=gmb;
-                else mb=fmb;
+
+                if(migration_style==SubelementMigrationEnum){
+                        if (phi<0.00000001) mb=gmb;
+                        else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
+                }
+                else if(migration_style==SubelementMigration2Enum){
+                        if(gllevelset>0.) mb=gmb;
+                        else mb=fmb;
+                }
+                else if(migration_style==SubelementMigration3Enum){
+                        if (phi<0.00000001) mb=gmb;
+                        else mb=gmb;
+                }
+                else{
+                        if(gllevelset>0) mb=gmb;
+                        else mb=fmb;
+                }
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
@@ -521 +550 @@
 
         /*Intermediaries */
-        int         migration_style,point1;
-        bool        mainlyfloating;
-        IssmDouble  Jdet,fraction1,fraction2;
-        IssmDouble  dt,gllevelset;
+        IssmDouble  Jdet,dt;
         IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.;
         IssmDouble* xyz_list = NULL;
-        Gauss*      gauss     = NULL;
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -545 +570 @@
         Input* thickness_input     = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
 
-        /*Recover portion of element that is grounded*/
-        phi=element->GetGroundedPortion(xyz_list);
-        if(migration_style==SubelementMigration2Enum){
-                gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
-                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
-           gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
-        }
-        else{
-                gauss = element->NewGauss(3);
-        }
-
         /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
@@ -566 +581 @@
                 gmb_input->GetInputValue(&gmb,gauss);
                 fmb_input->GetInputValue(&fmb,gauss);
-                gllevelset_input->GetInputValue(&gllevelset,gauss);
+                gllevelset_input->GetInputValue(&phi,gauss);
                 thickness_input->GetInputValue(&thickness,gauss);
-                if(migration_style==SubelementMigrationEnum){
-                        if (phi<0.00000001) mb=gmb;
-                        else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
-                }
-                else if(migration_style==SubelementMigration2Enum){
-                        if(gllevelset>0.) mb=gmb;
-                        else mb=fmb;
-                }
-                else if(migration_style==SubelementMigration3Enum){
-                        if (phi<0.00000001) mb=gmb;
-                        else mb=gmb;
-                }
-                else{
-                        if(gllevelset>0) mb=gmb;
-                        else mb=fmb;
-                }
+
+                if(phi>0.) mb=gmb;
+                else mb=fmb;
 
                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];