Changeset 13574

Timestamp:

10/09/12 21:05:46 (12 years ago)

Author:

Mathieu Morlighem

Message:

BUG: rewrote binary_search for double, there was a segfault in the previous implementation

Location:

issm/trunk-jpl/src/c/classes

Files:

• FemModel.cpp (modified) (1 diff)
• objects/Inputs/TransientInput.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -187 +187 @@
         /*Now delete: */
         delete profiler;
-
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp

@@ -429 +429 @@
         IssmDouble  alpha1,alpha2;
         
-        Input*  input=NULL;
-        Input*  input1=NULL;
-        Input*  input2=NULL;
+        Input *input  = NULL;
+        Input *input1 = NULL;
+        Input *input2 = NULL;
         
         bool    found=false;
@@ -440 +440 @@
          *fall within. Then interpolate the values on this interval: */
         found_offset=binary_search(&offset,intime,this->timesteps,this->numtimesteps);
-
-        if (found_offset==0){
-                if(intime<this->timesteps[0]){
-                        /*get values for the first time: */
-                        input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
-                        found=true;
-                }
-                else if(intime>this->timesteps[this->numtimesteps-1]){
-                        /*get values for the last time: */
-                        input=(Input*)((Input*)this->inputs->GetObjectByOffset(numtimesteps-1))->copy();
-                        found=true;
-                }
+        if (found_offset==0) _error_("Input not found (is TransientInput sorted ?)");
+
+        if (found_offset==-1){
+                /*get values for the first time: */
+                _assert_(intime<this->timesteps[0]);
+                input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy();
+                found=true;
+        }
+        else if(found_offset==3){
+                /*get values for the last time: */
+                _assert_(intime>=this->timesteps[this->numtimesteps-1]);
+                input=(Input*)((Input*)this->inputs->GetObjectByOffset(this->numtimesteps-1))->copy();
+                found=true;
+        }
+        else if(found_offset==2){
+                /*get values for this time: */
+                _assert_(intime==this->timesteps[offset]);
+                input=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy();
+                found=true;
         }
         else{
-                if (found_offset==1){
-                        /*intime is exactly equal to this->timesteps[offset]: */
-                        input=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy();
-                        found=true;
-                }
-                else if (found_offset==2){
-                        
-                        /*ok, we have the interval ]offset:offset+1[. Interpolate linearly for now: */
-                        deltat=this->timesteps[offset+1]-this->timesteps[offset];
-                        alpha2=(intime-this->timesteps[offset])/deltat;
-                        alpha1=(1.0-alpha2);
-
-                        input1=(Input*)this->inputs->GetObjectByOffset(offset); 
-                        input2=(Input*)this->inputs->GetObjectByOffset(offset+1);
-
-                        input=(Input*)input1->copy();
-                        input->Scale(alpha1);
-                        input->AXPY(input2,alpha2);
-
-                        found=true;
-                }
-                else _error_("binary_search returned the following value for found: " << found_offset<< "which is not supported yet!");
-        }
-
-        if(!found)_error_("did not find time interval on which to interpolate forcing values!");
+                /*get values between two times ]offset:offset+1[, Interpolate linearly*/
+                _assert_(intime>this->timesteps[offset] && intime<this->timesteps[offset+1]);
+                deltat=this->timesteps[offset+1]-this->timesteps[offset];
+                alpha2=(intime-this->timesteps[offset])/deltat;
+                alpha1=(1.0-alpha2);
+
+                input1=(Input*)this->inputs->GetObjectByOffset(offset); 
+                input2=(Input*)this->inputs->GetObjectByOffset(offset+1);
+
+                input=(Input*)input1->copy();
+                input->Scale(alpha1);
+                input->AXPY(input2,alpha2);
+
+                found=true;
+        }
 
         /*Assign output pointer*/