Changeset 17407

Timestamp:

03/11/14 07:22:01 (11 years ago)

Author:

jbondzio

Message:

CHG: extrapolate by diffusion or along normal by setting bool 'extrapolatebydiffusion'

File:

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

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

@@ -57 +57 @@
         if(VerboseSolution()) _printf0_("extrapolation of " << EnumToStringx(extvar_enum) << ": call computational core:\n");
         solutionsequence_linear(femmodel);
-        
-        save_results=true;
+
+        save_results=false;
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
@@ -79 +79 @@
         const int dim = 2;
         int        i,row,col,stabilization;
+        bool       extrapolatebydiffusion = true;
         IssmDouble Jdet,D_scalar,h;
         IssmDouble dlsf[dim],normal[dim];
@@ -90 +91 @@
         /*Initialize Element vector and other vectors*/
         ElementMatrix* Ke     = element->NewElementMatrix();
+        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
         IssmDouble     D[dim][dim];
 
         /*Retrieve all inputs and parameters*/
+        Input* lsf_slopex_input=element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
+        Input* lsf_slopey_input=element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
         element->GetVerticesCoordinates(&xyz_list);
+        h = element->CharacteristicLength();
 
         /* Start  looping on the number of gaussian points: */
@@ -102 +107 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                GetB(B,element,xyz_list,gauss);
                 GetBprime(Bprime,element,xyz_list,gauss);
-
+                
                 D_scalar=gauss->weight*Jdet;
 
-                for(row=0;row<dim;row++)
-                        for(col=0;col<dim;col++)
-                                D[row][col]=((row==col)?D_scalar:0.);
-
-                TripleMultiply(Bprime,dim,numnodes,1,
+                if(extrapolatebydiffusion){
+                        /* diffuse values outward */
+                        for(row=0;row<dim;row++)
+                                for(col=0;col<dim;col++)
+                                        if(row==col)
+                                                D[row][col] = D_scalar;
+                                        else
+                                                D[row][col] = 0.;
+
+                        TripleMultiply(Bprime,dim,numnodes,1,
                                         &D[0][0],dim,dim,0,
                                         Bprime,dim,numnodes,0,
                                         &Ke->values[0],1);
+                }
+                else{
+                        /* extrapolate values along normal */
+                        /* Get normal on ice boundary */
+                        lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
+                        lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
+                        norm_dlsf=0.;
+                        for(i=0;i<dim;i++)      norm_dlsf+=dlsf[i]*dlsf[i]; 
+                        norm_dlsf=sqrt(norm_dlsf); _assert_(norm_dlsf>0.);
+
+                        if(norm_dlsf>0.)
+                                for(i=0;i<dim;i++)      normal[i]=dlsf[i]/norm_dlsf;
+                        else
+                                for(i=0;i<dim;i++)      normal[i]=0.;
+                        
+                        for(row=0;row<dim;row++)
+                                for(col=0;col<dim;col++)
+                                        if(row==col)
+                                                D[row][col]=D_scalar*normal[row];
+                                        else
+                                                D[row][col]=0.;
+                        TripleMultiply(B,dim,numnodes,1,
+                                                &D[0][0],dim,dim,0,
+                                                Bprime,dim,numnodes,0,
+                                                &Ke->values[0],1);
+
+                        /* Stabilization *//*{{{*/
+                        stabilization=1;
+                        if (stabilization==0){/* no stabilization, do nothing*/}
+                        else if(stabilization==1){
+                                /* Artificial Diffusion */
+                                element->ElementSizes(&hx,&hy,&hz);
+                                h=sqrt( pow(hx*normal[0],2) + pow(hy*normal[1],2));
+                                kappa=h/2.+1.e-14; 
+                                D[0][0]=D_scalar*kappa;
+                                D[0][1]=0.;
+                                D[1][0]=0.;
+                                D[1][1]=D_scalar*kappa;
+                                TripleMultiply(Bprime,dim,numnodes,1,
+                                                        &D[0][0],dim,dim,0,
+                                                        Bprime,dim,numnodes,0,
+                                                        &Ke->values[0],1);
+                        }
+                        else if(stabilization==2){
+                                /*Streamline upwinding - do not use this for extrapolation: yields oscillating results due to smoothing along normal, not across */
+                                for(row=0;row<dim;row++)
+                                        for(col=0;col<dim;col++)
+                                                D[row][col]=h/(2.*1.)*normal[row]*normal[col];
+
+                                TripleMultiply(Bprime,dim,numnodes,1,
+                                                        &D[0][0],dim,dim,0,
+                                                        Bprime,dim,numnodes,0,
+                                                        &Ke->values[0],1);
+                        }/*}}}*/
+                }
         }/*}}}*/
 
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(Bprime);
         delete gauss;