Changeset 17274

Timestamp:

02/13/14 11:29:04 (11 years ago)

Author:

jbondzio

Message:

CHG: slope of levelset function for derivation of normal is now calculated once before extrapolation, since calculation in element can lead to normal = 0 for equal values on nodes.

Location:

issm/trunk-jpl/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/analyses/ExtrapolationAnalysis.cpp (modified) (3 diffs)
• c/analyses/L2ProjectionBaseAnalysis.cpp (modified) (2 diffs)
• c/classes/Elements/Penta.cpp (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (2 diffs)
• c/cores/cores.h (modified) (1 diff)
• c/cores/levelsetfunctionslope_core.cpp (added)
• c/cores/transient_core.cpp (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (1 diff)
• m/enum/EnumDefinitions.py (modified) (1 diff)
• m/enum/LevelsetfunctionSlopeXEnum.m (added)
• m/enum/LevelsetfunctionSlopeYEnum.m (added)

issm/trunk-jpl/src/c/Makefile.am ¶

@@ -421 +421 @@
                                         ./cores/meshdeformation_core.cpp\
                                         ./cores/damage_core.cpp\
+                                        ./cores/levelsetfunctionslope_core.cpp\
                                         ./analyses/DamageEvolutionAnalysis.cpp\
                                         ./solutionsequences/solutionsequence_damage_nonlinear.cpp\

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

@@ -74 +74 @@
         int        i,row,col,stabilization;
         IssmDouble Jdet,D_scalar,h;
-        IssmDouble dlevelset[dim],normal[dim];
-        IssmDouble norm_dlevelset;
+        IssmDouble dlsf[dim],normal[dim];
+        IssmDouble norm_dlsf;
         IssmDouble hx,hy,hz,kappa;
         IssmDouble* xyz_list = NULL;
@@ -89 +89 @@
 
         /*Retrieve all inputs and parameters*/
-        Input* levelset_input=element->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
+        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();
@@ -102 +103 @@
                 GetBprime(Bprime,element,xyz_list,gauss);
 
-                /* Get normal on node */
-                levelset_input->GetInputDerivativeValue(&dlevelset[0],xyz_list,gauss);
-                norm_dlevelset=0.;
-                for(i=0;i<dim;i++) norm_dlevelset+=dlevelset[i]*dlevelset[i]; 
-                norm_dlevelset=sqrt(norm_dlevelset)+1.e-14;
-                for(i=0;i<dim;i++) normal[i]=dlevelset[i]/norm_dlevelset;
-
+                /* 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;
 

issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp ¶

@@ -169 +169 @@
                 case BedSlopeXEnum:     input2 = basalelement->GetInput(BedEnum);     _assert_(input2); break;
                 case BedSlopeYEnum:     input2 = basalelement->GetInput(BedEnum);     _assert_(input2); break;
+                case LevelsetfunctionSlopeXEnum: input2 = basalelement->GetInput(MaskIceLevelsetEnum);     _assert_(input2); break;
+                case LevelsetfunctionSlopeYEnum: input2 = basalelement->GetInput(MaskIceLevelsetEnum);     _assert_(input2); break;
                 default: input = element->GetInput(input_enum);
         }
@@ -182 +184 @@
                 if(input2) input2->GetInputDerivativeValue(&slopes[0],xyz_list,gauss);
                 switch(input_enum){
-                        case SurfaceSlopeXEnum: case BedSlopeXEnum: value = slopes[0]; break;
-                        case SurfaceSlopeYEnum: case BedSlopeYEnum: value = slopes[1]; break;
+                        case SurfaceSlopeXEnum: case BedSlopeXEnum: case LevelsetfunctionSlopeXEnum: value = slopes[0]; break;
+                        case SurfaceSlopeYEnum: case BedSlopeYEnum: case LevelsetfunctionSlopeYEnum: value = slopes[1]; break;
                         default: input->GetInputValue(&value,gauss);
                 }

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶

@@ -1977 +1977 @@
                                 name==MaskIceLevelsetEnum ||
                                 name==IceMaskNodeActivationEnum ||
+                                name==LevelsetfunctionSlopeXEnum ||
+                                name==LevelsetfunctionSlopeYEnum ||
                                 name==GradientEnum ||
                                 name==OldGradientEnum  ||

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp ¶

@@ -910 +910 @@
         }
 
-        if(nrfrontnodes!=2){
-                _printf0_("in elt nr: " << this->Id() << "\n");
-                _printf0_("nrfrontnodes: " << nrfrontnodes << "\n");
-                for(i=0;i<NUMVERTICES;i++){
-                        _printf_("levelsetfunction[" << i << "]=" << levelset[i] << "; indicesfront: " << indicesfront[i] << ";\n");
-                }
-        }
         _assert_(nrfrontnodes==2);
 
@@ -1622 +1615 @@
                                 name==MaskIceLevelsetEnum ||
                                 name==IceMaskNodeActivationEnum ||
+                                name==LevelsetfunctionSlopeXEnum ||
+                                name==LevelsetfunctionSlopeYEnum ||
                                 name==SurfaceSlopeXEnum ||
                                 name==SurfaceSlopeYEnum ||

issm/trunk-jpl/src/c/cores/cores.h ¶

@@ -24 +24 @@
 void thermal_core(FemModel* femmodel);
 void surfaceslope_core(FemModel* femmodel);
+void levelsetfunctionslope_core(FemModel* femmodel);
 void bedslope_core(FemModel* femmodel);
 void meshdeformation_core(FemModel* femmodel);

issm/trunk-jpl/src/c/cores/transient_core.cpp ¶

@@ -125 +125 @@
                 if(islevelset){
                         if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
+                        
+                        /* get slope of lsf for computation of normal on ice domain*/
+                        levelsetfunctionslope_core(femmodel);
 
-                        /* extrapolate required variables */
+                        /* extrapolate along normal */
                         Analysis* extanalysis = new ExtrapolationAnalysis();
                         const int nvars=2;

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -679 +679 @@
         ExtrapolationVariableEnum,
         IceMaskNodeActivationEnum,
+        LevelsetfunctionSlopeXEnum,
+        LevelsetfunctionSlopeYEnum,
         /*}}}*/
         MaximumNumberOfDefinitionsEnum

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -638 +638 @@
                 case ExtrapolationVariableEnum : return "ExtrapolationVariable";
                 case IceMaskNodeActivationEnum : return "IceMaskNodeActivation";
+                case LevelsetfunctionSlopeXEnum : return "LevelsetfunctionSlopeX";
+                case LevelsetfunctionSlopeYEnum : return "LevelsetfunctionSlopeY";
                 case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
                 default : return "unknown";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -653 +653 @@
               else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
               else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
+              else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
+              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
               else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
          else stage=7;

issm/trunk-jpl/src/m/enum/EnumDefinitions.py ¶

@@ -630 +630 @@
 def ExtrapolationVariableEnum(): return StringToEnum("ExtrapolationVariable")[0]
 def IceMaskNodeActivationEnum(): return StringToEnum("IceMaskNodeActivation")[0]
+def LevelsetfunctionSlopeXEnum(): return StringToEnum("LevelsetfunctionSlopeX")[0]
+def LevelsetfunctionSlopeYEnum(): return StringToEnum("LevelsetfunctionSlopeY")[0]
 def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]