-
Notifications
You must be signed in to change notification settings - Fork 3
/
TransientArrayParam.cpp
147 lines (125 loc) · 4.33 KB
/
TransientArrayParam.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
/*!\file TransientArrayParam.c
* \brief: implementation of the TransientArrayParam object
*/
/*header files: */
/*{{{*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
#include "../classes.h"
#include "../../shared/shared.h"
/*}}}*/
/*TransientArrayParam constructors and destructor*/
TransientArrayParam::TransientArrayParam(){/*{{{*/
return;
}
/*}}}*/
TransientArrayParam::TransientArrayParam(int in_enum_type,IssmDouble* in_values,IssmDouble* in_time,bool interpolation_on,bool cycle_in,int in_N,int in_M){/*{{{*/
_assert_(in_values && in_time);
this->enum_type=in_enum_type;
this->M=in_M; //Number of rows
this->N=in_N; //Number of timesteps
this->interpolation=interpolation_on;
this->cycle=cycle_in;
this->values=xNew<IssmDouble>(M*N);
xMemCpy<IssmDouble>(this->values,in_values,M*N);
this->timesteps=xNew<IssmDouble>(N);
xMemCpy<IssmDouble>(this->timesteps,in_time,N);
}
/*}}}*/
TransientArrayParam::~TransientArrayParam(){/*{{{*/
xDelete<IssmDouble>(values);
xDelete<IssmDouble>(timesteps);
}/*}}}*/
/*Object virtual functions definitions:*/
Param* TransientArrayParam::copy() {/*{{{*/
return new TransientArrayParam(this->enum_type,this->values,this->timesteps,this->interpolation,this->cycle,this->M,this->N);
}
/*}}}*/
void TransientArrayParam::DeepEcho(void){/*{{{*/
_printf_("TransientArrayParam:\n");
_printf_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
_printf_(" number of time steps: " << this->N << "\n");
_printf_(" number of rows: " << this->M << "\n");
for(int i=0;i<this->N;i++){
_printf_(" time: " << this->timesteps[i] << "\n");
for(int k=0;k<this->M;k++){
_printf_(" values: " << this->values[k*N+i] << "\n");
}
_printf_("\n");
}
}
/*}}}*/
void TransientArrayParam::Echo(void){/*{{{*/
_printf_("TransientArrayParam:\n");
_printf_(" enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
_printf_(" size: " << this->M << " by " << this->N << "\n");
}
/*}}}*/
int TransientArrayParam::Id(void){ return -1; }/*{{{*/
/*}}}*/
void TransientArrayParam::Marshall(MarshallHandle* marshallhandle){ /*{{{*/
int object_enum = TransientArrayParamEnum;
marshallhandle->call(object_enum);
marshallhandle->call(this->enum_type);
marshallhandle->call(this->interpolation);
marshallhandle->call(this->cycle);
marshallhandle->call(this->M);
marshallhandle->call(this->N);
if(marshallhandle->OperationNumber()==MARSHALLING_LOAD){
values = xNew<IssmDouble>(M*N);
timesteps = xNew<IssmDouble>(N);
}
marshallhandle->call(this->values,M*N);
marshallhandle->call(this->timesteps,N);
}/*}}}*/
int TransientArrayParam::ObjectEnum(void){/*{{{*/
return TransientArrayParamEnum;
}/*}}}*/
/*TransientArrayParam virtual functions definitions: */
void TransientArrayParam::GetParameterValue(IssmDouble* pdouble,int row,IssmDouble time){/*{{{*/
IssmDouble output;
bool found;
_assert_(row>=0 && row<this->M);
if(this->cycle) _error_("not implemented yet");
/*Ok, we have the time and row, go through the timesteps, and figure out which interval we
*fall within. Then interpolate the values on this interval: */
if(time<this->timesteps[0]){
/*get values for the first time: */
output=this->values[row*this->N];
found=true;
}
else if(time>this->timesteps[this->N-1]){
/*get values for the last time: */
output=this->values[(row+1)*this->N-1];
found=true;
}
else{
/*Find which interval we fall within: */
for(int i=0;i<this->N;i++){
if(time==this->timesteps[i]){
/*We are right on one step time: */
output = this->values[row*this->N+i];
found=true;
break; //we are done with the time interpolation.
}
else{
if(this->timesteps[i]<time && time<this->timesteps[i+1]){
/*ok, we have the interval [i:i+1]. Interpolate linearly for now: */
IssmDouble deltat = this->timesteps[i+1]-this->timesteps[i];
IssmDouble alpha = (time-this->timesteps[i])/deltat;
if(this->interpolation==true) output=(1.0-alpha)*this->values[row*this->N+i] + alpha*this->values[row*this->N+i+1];
else output=this->values[row*this->N+i];
found=true;
break;
}
else continue; //keep looking on the next interval
}
}
}
if(!found)_error_("did not find time interval on which to interpolate values");
*pdouble=output;
}
/*}}}*/