Changeset 13577
- Timestamp:
- 10/10/12 07:39:19 (12 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/objects/Inputs/TransientInput.cpp
r13574 r13577 426 426 Input* TransientInput::GetTimeInput(IssmDouble intime){ 427 427 428 IssmDouble deltat; 429 IssmDouble alpha1,alpha2; 428 IssmDouble deltat; 429 IssmDouble alpha1,alpha2; 430 int found; 431 int offset; 430 432 431 433 Input *input = NULL; … … 433 435 Input *input2 = NULL; 434 436 435 bool found=false; 436 int found_offset=0; 437 int offset; 438 439 /*Ok, we have the time, go through the timesteps, and figure out which interval we 437 /*go through the timesteps, and figure out which interval we 440 438 *fall within. Then interpolate the values on this interval: */ 441 found _offset=binary_search(&offset,intime,this->timesteps,this->numtimesteps);442 if (found_offset==0) _error_("Input not found (is TransientInput sorted ?)");443 444 if ( found_offset==-1){439 found=binary_search(&offset,intime,this->timesteps,this->numtimesteps); 440 if(!found) _error_("Input not found (is TransientInput sorted ?)"); 441 442 if (offset==-1){ 445 443 /*get values for the first time: */ 446 444 _assert_(intime<this->timesteps[0]); 447 445 input=(Input*)((Input*)this->inputs->GetObjectByOffset(0))->copy(); 448 found=true; 449 } 450 else if(found_offset==3){ 446 } 447 else if(offset==(this->numtimesteps-1)){ 451 448 /*get values for the last time: */ 452 _assert_(intime>=this->timesteps[this->numtimesteps-1]); 453 input=(Input*)((Input*)this->inputs->GetObjectByOffset(this->numtimesteps-1))->copy(); 454 found=true; 455 } 456 else if(found_offset==2){ 457 /*get values for this time: */ 458 _assert_(intime==this->timesteps[offset]); 449 _assert_(intime>=this->timesteps[offset]); 459 450 input=(Input*)((Input*)this->inputs->GetObjectByOffset(offset))->copy(); 460 found=true;461 451 } 462 452 else{ 463 /*get values between two times ]offset:offset+1[, Interpolate linearly*/464 _assert_(intime> this->timesteps[offset] && intime<this->timesteps[offset+1]);453 /*get values between two times [offset:offset+1[, Interpolate linearly*/ 454 _assert_(intime>=this->timesteps[offset] && intime<this->timesteps[offset+1]); 465 455 deltat=this->timesteps[offset+1]-this->timesteps[offset]; 466 456 alpha2=(intime-this->timesteps[offset])/deltat; … … 473 463 input->Scale(alpha1); 474 464 input->AXPY(input2,alpha2); 475 476 found=true;477 465 } 478 466 -
issm/trunk-jpl/src/c/shared/Sorting/binary_search.cpp
r13575 r13577 64 64 return found; 65 65 } /*}}}*/ 66 int binary_search(int* poffset,double target,double* list,int num_doubles){ /*{{{*/ 66 int binary_search(int* poffset,double target,double* list,int length){ /*{{{*/ 67 /* 68 * l[0] l[1] l[2] l[n] l[n+1] l[length-1] 69 * <-------+-----+-----+-----+ ... +-----+........+--------------> 70 * offset: -1 0 1 2 n length -1 71 * 72 * offset = -1 target < list[0] 73 * offset = n list[n] <= target < list[n+1] 74 * offset = length-1 list[length-1] <= target 75 */ 67 76 68 77 /*output: */ 69 78 int offset = 0; 70 int found = 0; /*found = 0: not found. 71 found = -1: found, and target is < first element 72 found = 1: found, and target is == to value at offset 73 found = 2: found, and target is > to value at offset and < to value at offset+1 74 found = 3: found, and target is >= last value */ 79 int found = 0; 75 80 76 81 /*intermediary: */ 77 82 int n0 = 0; 78 int n1 = int( num_doubles/2);79 int n2 = num_doubles-1;83 int n1 = int(length/2); 84 int n2 = length-1; 80 85 81 86 if(target<list[n0]){ 82 found = -1;87 found = 1; 83 88 offset = -1; 84 89 } 85 90 else if(target>=list[n2]){ 86 found = 3;87 offset = n2-1;91 found = 1; 92 offset = length-1; 88 93 } 89 94 else{ 90 while(!found){95 for(;;){ 91 96 /*did we find the target?*/ 92 97 if(list[n1]<=target && list[n1+1]>target){ 93 found = 1;98 found = 1; 94 99 offset = n1; 95 100 break; 96 101 } 97 if(target < list[n1]){102 else if(target < list[n1]){ 98 103 n2 = n1; 99 104 n1 = n0 + int((n2-n0)/2); … … 104 109 } 105 110 } 106 107 //did we find an exact target?108 if(list[n1]==target) found = 2;109 111 } 110 112 111 /*Assign output pointer s:*/113 /*Assign output pointer and return*/ 112 114 *poffset=offset; 113 114 /*Return result: */115 115 return found; 116 116 } /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.