Permanantly (hopefully) fix out-of-bounds errors for boundary
authorOrion Lawlor <olawlor@acm.org>
Thu, 13 Jan 2005 00:50:07 +0000 (00:50 +0000)
committerOrion Lawlor <olawlor@acm.org>
Thu, 13 Jan 2005 00:50:07 +0000 (00:50 +0000)
condition maintainance.  If this doesn't work, just hardcode
prop=1.0 and slope=0.0.

examples/fem/crack2D/crack.h
examples/fem/crack2D/node.C

index d31e18c435cc696dfd0d5299f4be3c660bb276ae..ef7aa607611c6d08b9b1bfe3c420d6f4a76c9ab2 100644 (file)
@@ -98,9 +98,16 @@ struct ConfigurationData {
   double delta;            //timestep
   double delta2;           //delta squared times 2????
   int numProp;             //number of proportionality constants
-  int *ts_proportion;      //time step to apply proportionality constant at
+  int *ts_proportion;      //1-based time step to apply proportionality constant at
   double *proportion;      //boundary load proportionality constant
 
+  /**
+   Return the fraction of the boundary conditions (prop) and time-rate of
+   change of boundary conditions (slope) to apply during this 0-based timestep.  
+   Applies linear interpolation to the ts_proportion and proportion arrays, above.
+  */
+  void getBoundaryConditionScale(int timestep,double *prop,double *slope) const;
+
   /// Material formulation
   int lin;                 //elastic formulation, 0=nl, 1=linear
   int ncoh;                //type of cohesive law, 0=exp., 1=linear
@@ -153,7 +160,6 @@ extern "C" void pupMesh(pup_er p,MeshData *mesh); //For migration
 
 //Node physics: in node.C
 struct NodeSlope {
-  int kk; //Index into config.ts_proportion array
   double prop; //Proportion of boundary conditions to apply
   double slope; //Slope of prop
 };
index 580a9d70fb0f87e9c447924bc8d46afe8073606f..c55426314e8619956fa9fb3ad376f6f9e79ce89b 100644 (file)
@@ -6,7 +6,6 @@
 #include "crack.h"
 
 void nodeSetup(NodeSlope *sl) {
-  sl->kk=-1;
   sl->slope=sl->prop=0;
 }
 
@@ -21,24 +20,38 @@ void nodeBeginStep(MeshData *mesh)
   }
 }
 
+/**
+ Return the fraction of the boundary conditions (prop) and time-rate of
+ change of boundary conditions (slope) to apply during this 0-based timestep.  
+ Applies linear interpolation to the ts_proportion and proportion arrays, above.
+*/
+void ConfigurationData::getBoundaryConditionScale(int timestep,double *prop,double *slope) const
+{
+    timestep--; /* to 0-based */
+    /* Clamp out-of-bounds timesteps */
+    if (timestep<=ts_proportion[0]) {
+      *prop=proportion[0]; *slope=0; return;
+    }
+    if (timestep>=ts_proportion[numProp-1]) {
+      *prop=proportion[numProp-1]; *slope=0; return;
+    }
+    /* Otherwise linearly interpolate between proportion "cur" and "cur+1" */
+    int cur=0; /* index into ts_proportion and proportion arrays */
+    while (timestep>=ts_proportion[cur+1]) cur++;
+    /* assert: 0<=cur && cur<numProp 
+      and  ts_proportion[cur]<=timestep && timestep<ts_proportion[cur+1] */
+    
+    *prop = proportion[cur];
+    double timeSlope=1.0/(ts_proportion[cur+1]-ts_proportion[cur]);
+    *slope = 1/delta * (proportion[cur+1]-proportion[cur])*timeSlope;
+    *prop = proportion[cur]+(double)(timestep-ts_proportion[cur])*(*slope)*delta;
+}
+
+
 void nodeFinishStep(MeshData *mesh, NodeSlope *sl,int tstep)
 {
-  tstep++; /* subtle: timesteps are 1-based in the .inp file, so make tstep 1-based */
-  // Slowly ramp up boundary conditions:
-  if (config.ts_proportion[sl->kk+1] == tstep)
-  {
-      sl->kk++;
-      sl->prop = config.proportion[sl->kk];
-      sl->slope = (config.proportion[sl->kk+1]-sl->prop)/config.delta;
-      sl->slope /= (double) (config.ts_proportion[sl->kk+1]- config.ts_proportion[sl->kk]);
-  }
-  else
-  {
-      sl->prop = (double)(tstep - config.ts_proportion[sl->kk])*
-                    sl->slope*config.delta+config.proportion[sl->kk];
-  }
-  double slope=sl->slope;
-  double prop=sl->prop;
+  double prop,slope;
+  config.getBoundaryConditionScale(tstep,&prop,&slope);
 
   // Update each node:
   for(int idx=0; idx<mesh->nn; idx++) {