Changeset 21759 for issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/ExtrapolationAnalysis.cpp

Timestamp:

06/07/17 10:50:54 (8 years ago)

Author:

Eric.Larour

Message:

CHG: merged branch back to trunk-jpl 21754.

File:

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

issm/branches/trunk-larour-NatGeoScience2016/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;
+        bool          extrapolatebydiffusion = true;
+        int           dim, domaintype, extrapolationcase;
+        int           i,row,col,stabilization;
+        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*/
@@ -126 +128 @@
         workelement->GetVerticesCoordinates(&xyz_list);
 
+        /* In 3d we are going to extrude using horizontal diffusion. Since the layers are 3d,
+         * we change the geometry of the element to make it flat. That way, the diffusion is
+         * handled along the layers
+         */
+        if(element->ObjectEnum()==PentaEnum){
+                for(i=0;i<3;i++) xyz_list[3*i+2] = 0.;
+                for(i=3;i<6;i++) xyz_list[3*i+2] = 1.;
+        }
+
         /* 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);
 
@@ -137 +148 @@
                 D_scalar=gauss->weight*Jdet;
 
-                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 +186 @@
                                                 &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 +206 @@
                                                         &Ke->values[0],1);
                         }
-                        /*}}}*/
                 }
-        }/*}}}*/
+        }
 
         /*Clean up and return*/
@@ -217 +222 @@
 }/*}}}*/
 ElementVector* ExtrapolationAnalysis::CreatePVector(Element* element){/*{{{*/
-
-        /*Intermediaries */
-        Element* workelement=NULL;
-
-        /*Get problem dimension*/
-        int extrapolationcase=GetExtrapolationCase(element);
-        switch(extrapolationcase){
-                case 0: 
-                        if(!element->IsOnBase()) return NULL; 
-                        workelement = element->SpawnBasalElement(); 
-                        break;
-                case 1: case 2: case 3: workelement=element; break;
-        }
-
-        /*Fetch number of nodes */
-        int numnodes = workelement->GetNumberOfNodes();
-
-        /*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;};
-        return pe;
+        return NULL;
+
 }/*}}}*/
 void           ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss, int dim){/*{{{*/
@@ -326 +309 @@
         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;
         }