Changeset 21401

Timestamp:

11/18/16 13:37:29 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: cosmetics

File:

• issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (modified) (9 diffs)

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

@@ -84 +84 @@
 
         /*Intermediaries */
-        int     dim, domaintype, extrapolationcase;
-        int     i,row,col,stabilization;
-        bool    extrapolatebydiffusion;
-        IssmDouble Jdet,D_scalar,h;
-        IssmDouble norm_dlsf;
-        IssmDouble hx,hy,hz,kappa;
-        IssmDouble*     xyz_list = NULL;
-        Element*                workelement=NULL;
+        int           dim, domaintype, extrapolationcase;
+        int           i,row,col,stabilization;
+        bool          extrapolatebydiffusion;
+        IssmDouble  Jdet,D_scalar,h;
+        IssmDouble  norm_dlsf;
+        IssmDouble  hx,hy,hz,kappa;
+        IssmDouble *xyz_list    = NULL;
+        Element    *workelement = NULL;
 
         /*Get problem case*/
@@ -102 +102 @@
                 case 1: case 2: case 3: workelement=element; break;
         }
+
         /* get extrapolation dimension */
         workelement->FindParam(&domaintype,DomainTypeEnum);
         switch(domaintype){
-                case Domain2DverticalEnum: dim=1; break;
+                case Domain2DverticalEnum:   dim=1; break;
                 case Domain2DhorizontalEnum: dim=2; break;
-                case Domain3DEnum: dim=3; break;
+                case Domain3DEnum:           dim=3; break;
+      default: _error_("not supported yet");
         }
 
@@ -117 +119 @@
         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*             D         = xNew<IssmDouble>(dim*dim);
-        IssmDouble*             dlsf  = xNew<IssmDouble>(dim);
-        IssmDouble*             normal= xNew<IssmDouble>(dim);
+        IssmDouble*             dlsf   = xNew<IssmDouble>(dim);
+        IssmDouble*             normal = xNew<IssmDouble>(dim);
+   IssmDouble*          D           = xNewZeroInit<IssmDouble>(dim*dim);
 
         /*Retrieve all inputs and parameters*/
@@ -128 +130 @@
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=workelement->NewGauss(2);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){/*{{{*/
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
                 gauss->GaussPoint(ig);
 
@@ -139 +141 @@
                 extrapolatebydiffusion=true;
                 if(extrapolatebydiffusion){
-                        /* diffuse values outward */
-                        for(row=0;row<dim;row++)
-                                for(col=0;col<dim;col++)
-                                        if(row==col && row<2) //extrapolate only in xy-plane
-                                                D[row*dim+col] = D_scalar;
-                                        else
-                                                D[row*dim+col] = 0.;
+
+                        /* diffuse values outward only along the xy-plane*/
+         for(int i=0;i<2;i++) D[i*dim+i] = D_scalar;
 
                         TripleMultiply(Bprime,dim,numnodes,1,
@@ -180 +178 @@
                                                 &Ke->values[0],1);
 
-                        /* stabilization *//*{{{*/
+                        /* stabilization */
                         /* do not use streamline upwinding for extrapolation: it yields oscillating results due to diffusion along normal direction, but none across */
                         stabilization=1;
@@ -200 +198 @@
                                                         &Ke->values[0],1);
                         }
-                        /*}}}*/
                 }
-        }/*}}}*/
+        }
 
         /*Clean up and return*/
@@ -236 +233 @@
         /*Initialize Element vector*/
         ElementVector* pe = workelement->NewElementVector();
-        for(int i=0;i<numnodes;i++) 
-                pe->values[i]=0.; 
-
-        if(extrapolationcase==0){workelement->DeleteMaterials(); delete workelement;};
+
+        if(extrapolationcase==0){
+      workelement->DeleteMaterials();
+      delete workelement;
+   }
+
         return pe;
 }/*}}}*/
@@ -326 +325 @@
         int domaintype, extrapolationvariable;
         int extrapolationcase;
+
         element->FindParam(&domaintype,DomainTypeEnum);
         switch(domaintype){
-                case Domain2DverticalEnum: extrapolationcase=0; break;
-                case Domain2DhorizontalEnum: extrapolationcase=1;break;
-                case Domain3DEnum:  
+                case Domain2DverticalEnum:   extrapolationcase=0; break;
+                case Domain2DhorizontalEnum: extrapolationcase=1; break;
+                case Domain3DEnum: 
                         element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
-                        if(extrapolationvariable==ThicknessEnum) extrapolationcase=2; // scalar fields that are constant along z-axis
-                        else extrapolationcase=3; // scalar fields that vary along z-axis
+                        if(extrapolationvariable==ThicknessEnum){
+            extrapolationcase=2; // scalar fields that are constant along z-axis
+         }
+                        else{
+            extrapolationcase=3; // scalar fields that vary along z-axis
+         }
                         break;
         }