Changeset 17429 for issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp

Timestamp:

03/13/14 15:32:02 (11 years ago)

Author:

jbondzio

Message:

BUG: deactivate nodes on non-basal layer for 2d solving process of extrapolation and LSM

File:

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

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

@@ -29 +29 @@
                 }
         }
-        
         if(iomodel->meshtype==Mesh3DEnum){
                 iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
+                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
         }
 }
@@ -86 +86 @@
 
         /*Intermediaries */
-        int               meshtype,dim;
+        int                meshtype,dim;
         int        i,row,col,stabilization;
         bool       extrapolatebydiffusion = true;
@@ -110 +110 @@
         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
-        IssmDouble*             dlsf   = xNew<IssmDouble>(dim);
-        IssmDouble*             normal = xNew<IssmDouble>(dim);
-        IssmDouble*    D                 = xNew<IssmDouble>(dim*dim);
+        IssmDouble*             D         = xNew<IssmDouble>(dim*dim);
+        IssmDouble*             dlsf  = xNew<IssmDouble>(dim);
+        IssmDouble*             normal= xNew<IssmDouble>(dim);
 
         /*Retrieve all inputs and parameters*/
@@ -162 +162 @@
                                 for(col=0;col<dim;col++)
                                         if(row==col)
-                                                D[row*dim+col] = D_scalar*normal[row];
+                                                D[row*dim+col]=D_scalar*normal[row];
                                         else
-                                                D[row*dim+col] = 0.;
+                                                D[row*dim+col]=0.;
                         TripleMultiply(B,dim,numnodes,1,
                                                 D,dim,dim,0,
@@ -181 +181 @@
                                         for(col=0;col<dim;col++)
                                                 if(row==col)
-                                                        D[row*dim+col] = D_scalar*kappa;
+                                                        D[row*dim+col]=D_scalar*kappa;
                                                 else
-                                                        D[row*dim+col] = 0.;
+                                                        D[row*dim+col]=0.;
                                 TripleMultiply(Bprime,dim,numnodes,1,
                                                         D,dim,dim,0,
@@ -207 +207 @@
         xDelete<IssmDouble>(B);
         xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(D);
         xDelete<IssmDouble>(dlsf);
         xDelete<IssmDouble>(normal);
+        delete gauss;
         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-        delete gauss;
         return Ke;
 
@@ -224 +225 @@
         /*Fetch number of nodes */
         int numnodes = basalelement->GetNumberOfNodes();
-        basalelement->FindParam(&meshtype,MeshTypeEnum);
 
         /*Initialize Element vector*/
@@ -230 +230 @@
         for(i=0;i<numnodes;i++) 
                 pe->values[i]=0.; 
+
+        basalelement->FindParam(&meshtype,MeshTypeEnum);
         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         return pe;
@@ -326 +328 @@
                 node=element->GetNode(in);
                 levelset_input->GetInputValue(&phi,gauss);
-                if(phi<=0.){
-                        /* if ice, set dirichlet BC */
-                        extvar_input->GetInputValue(&value,gauss);
-                        node->ApplyConstraint(1,value);
-                }
-                else {
-                        /* no ice, set no spc */
-                        node->DofInFSet(0); 
+                if(node->IsActive()){
+                        if(phi<=0.){
+                                /* if ice, set dirichlet BC */
+                                extvar_input->GetInputValue(&value,gauss);
+                                node->ApplyConstraint(1,value);
+                        }
+                        else {
+                                /* no ice, set no spc */
+                                node->DofInFSet(0); 
+                        }
                 }
         }