Changeset 17363

Timestamp:

02/27/14 13:16:16 (11 years ago)

Author:

jbondzio

Message:

CHG: Instead of extrapolation along normal, diffuse variable outward. Once initialization of levelset function to signed distance function is working, revert this step.

File:

• issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (modified) (3 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;
+        if(save_results){
+                if(VerboseSolution()) _printf0_("   saving results\n");
+                int outputs[2] = {VxEnum,VyEnum};
+                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+        }
 
 }/*}}}*/
@@ -84 +90 @@
         /*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: */
@@ -100 +102 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                GetB(B,element,xyz_list,gauss);
                 GetBprime(Bprime,element,xyz_list,gauss);
 
-                /* 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.);
-                for(i=0;i<dim;i++)
-                        normal[i]=dlsf[i]/norm_dlsf;
-                
                 D_scalar=gauss->weight*Jdet;
 
                 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[row][col]=((row==col)?D_scalar:0.);
+
+                TripleMultiply(Bprime,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.; 
-                        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;